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ABSTRACT 

A theoretical model framework of spherical symmetry is presented for a composite 
astrophysical system of two polytropic fluids coupled together by gravity to explore 
large-scale shocks and flow dynamics in clusters of galaxies or in globular clusters. The 
existence of such large-scale shocks in clusters of galaxies as inferred by high-resolution 
X-ray and radio imaging observations implies large-scale systematic flows that are 
beyond usual static models for clusters of galaxies. Here, we explore self-similar two- 
fluid flow solutions with shocks for a hot polytropic gas flow in a cluster of galaxies in 
the presence of a massive dark matter (DM) flow after the initiation of a gravitational 
core collapse or a central AGN activity or a large-scale merging process. In particular, 
the possibility of DM shocks or sharp jumps of mass density and of velocity dispersion 
in dark matter halo is discussed and such DM shocks might be detectable through 
gravitational lensing effects. To examine various plausible scenarios for clusters of 
galaxies, we describe three possible classes of shock flows within our model framework 
for different types of temperature, density and flow speed profiles. Depending upon 
sensible model parameters and shock locations, the hot ICM and DM halo may have 
various combinations of asymptotic behaviours of outflow, breeze, inflow, contraction 
or static envelopes at large radii at a given time. We refer to asymptotic outflows of 
hot ICM at large radii as the galaxy cluster wind. As a result of such galaxy cluster 
winds and simultaneous contractions of DM halo during the course of galaxy cluster 
evolution, there would be less hot ICM within clusters of galaxies as compared to 
the average baryon fraction in the Universe. Physically, it is then expected that such 
'missing baryons' with lower temperatures reside in the periphery of galaxy clusters 
on much larger scales. Based on our model analysis, we also predict a limiting (the 
steepest) radial scaling form for mass density profiles of r~ 3 within clusters of galaxies. 

Key words: dark matter — galaxies: clusters: general — gravitation — hydrody- 
namics — shock waves — X-rays: galaxies: clusters 



1 INTRODUCTION 

Extensive X-ray observations have revealed that almost 
completely ionized hot gas medium permeates within clus- 
ters of galaxies with a typical temperature range of ~ 
10 7 - 10 s K and a range of typical electron number den- 
sity ~ 10~ 2 — 10 -4 cm -3 (e.g., Cavaliere & Fusco-Femiano 
1978; Sarazin 1988; Fabian 1994). Clusters of galaxies are 
largely gravitationally bound systems on spatial scales of 
several Mpcs; together with the strong evidence of high ve- 
locity dispersions of galaxies (e.g., ~ 700 — 1000 km s _1 ), 
hot X-ray emitting thermal electron gas, and gravitational 
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lensing effects, we have realized the presence of massive dark 
matter halo within clusters of galaxies. Physical properties 
of galaxy clusters with static models of spherical symmetry 
have been extensively studied in the past (e.g., Lea 1975; 
Cavaliere & Fusco-Femiano 1976, 1978; Sarazin & Bahcall 
1977; Sarazin 1988; Fabian 1994; Carilli & Taylor 2002; Voit 
2005). 

In the past several years, high-resolution X-ray imaging 
observations have revealed density, temperature and pres- 
sure jumps in the hot intracluster medium (ICM) within 
clusters of galaxies (e.g., Fabian et al. 2003; Nulsen et al. 
2005a; Nulsen et al. 2005b; McNamara et al. 2005), indicat- 
ing that these structures are likely large-scale shocks rather 
than cold fronts (e.g., Sanders & Fabian 2006). Active galac- 
tic nuclei (AGNs) (e.g., Nulsen et al. 2005a; Nulsen et al. 
2005b; McNamara et al. 2005) and merging galaxies (e.g., 
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Markevitch & Vikhlinin 2001; Markevitch et al. 2002; Gabici 
& Blasi 2003; Markevitch et al. 2005) are proposed to be the 
driving force and energy source of these large-scale shocks. 
In addition, large-scale sound waves in clusters of galaxies 
have been proposed by Sanders & Fabian (2007) to explain 
the observed quasi-concentric ripples in surface brightness 
of X-ray emissions. Once shocks are identified in clusters of 
galaxies, there must be large-scale flows involved. In other 
words, these clusters of galaxies cannot be really static on 
large scales at least in the spatial region where shocks are 
presumably identified. Relative to the cluster centre, radial 
distances of shocks observed vary from tens of kpcs (e.g., 
McNamara et al. 2005; Nulsen et al. 2005b) to several Mpcs 
(e.g., A3667 in Rotteringer et al. 1997; A3376 in Bagchi et al. 
2006; galaxy clusters A786, A2255, A2256 in Ensslin et al. 
1998); at yet smaller radii around the centre, there could be 
emerging shocks that may not be easily identified. We take 
the point of view that these shocks are moving in hot ICM 
and we simply catch them at different epochs of evolution. 
These shocks may all be born somewhere around the central 
region and travel outwards to the locations we observe at the 
present epoch. For those clusters of galaxies without shock 
signatures, one possibility is that shocks have occurred in 
the distant past and have disappeared after their energies 
were dissipated during the process of propagation. So large- 
scale shocks and flows may well be common phenomena in 
clusters of galaxies. 

On much smaller scales compared to those of clusters 
of galaxies, X-rays have been also observed in globular clus- 
ters (e.g., Verbunt et al. 1984) and these emissions are in- 
terpreted by some as associated with flowing gas towards a 
black hole residing in the centre of globular clusters (e.g., 
Silk & Arons 1975). Moreover, outflows of gas materials 
from globular clusters have also been discussed in the litera- 
ture (e.g., VandenBerg 1978). Fully relaxed globular clusters 
can be well treated as spherically symmetric (e.g., Harris & 
Racine 1979). One may view the collection of stars as one 
'fluid' and the tenuous gas as another fluid; these two flu- 
ids are coupled together by gravity on large scales. Here, we 
note globular clusters in passing and will mainly focus on 
large-scale self-similar dynamics for clusters of galaxies. 

With these two classes of astrophysical systems in mind, 
we develop a theoretical model framework of spherical sym- 
metry to study dynamic behaviours of hot ICM and dark 
matter halo using the two- fluid approximation (e.g., Lou 
2005), where the two polytropic fluids are coupled together 
by gravity. We should note that here the notion of a poly- 
tropic fluid is fairly general in the sense of specific entropy 
conservation along streamlines (Lou & Cao 2007). 

In the context of large-scale structure formation, ex- 
tensive numerical works have been carried out to simulate 
the formation of galaxy clusters in the expanding universe, 
providing information for the hot ICM and dark matter 
halo (see, e.g., Bertschinger 1998 for a review on numerical 
simulations of structure formation in the universe). Evrard 
(1990) and Thomas & Couchman (1992) simulated prop- 
erties (such as the number density and temperature pro- 
files) of a hot gas in the presence of a dark matter halo. 
In Katz & White (1993) and Frenk et al. (1996), the ra- 
dial cooling process is simulated. In particular, Evrard et al. 
(1994) simulated the formation of galaxies with two gravi- 
tationally coupled fluids representing dark matter halo and 



baryon matter, which is similar in essence to the approxi- 
mation adopted in our semi-analytical model for clusters of 
galaxies but on much larger scales. 

Various self-similar solutions describing hydrodynamic 
processes of a single self-gravitational isothermal or poly- 
tropic gas under spherical symmetry have been investigated 
previously in contexts of star formation (e.g., Larson 1969a; 
Larson 1969b; Shu 1977; Hunter 1977; Shu et al. 1987; Lou 
& Shen 2004). Very recently, asymptotic behaviours of novel 
quasi-static solutions in a single polytropic gas sphere with 
self-gravity have been reported by Lou & Wang (2006, 2007) 
and was utilized to model rebound (MHD) shocks in super- 
novae (Wang & Lou 2007). For an astrophysical system of 
two fluids coupled by gravity, we can systematically extend 
these self-similar solutions, especially the new quasi-static 
solution, which may be used to describe behaviours of hot 
ICM and dark matter halo in clusters of galaxies. Except 
for the gravitational effect in the Newtonian sense, nothing 
else is known about dark matter particles at present. Us- 
ing the coupled two-fluid model, we might be able to learn 
physical properties of dark matter halo through detectable 
diagnostics of hot ICM and of gravitational lensing effects. 

For clusters of galaxies, there is an outstanding prob- 
lem of 'missing baryons'. Extensive X-ray observations have 
indicated that the baryon mass fraction in clusters of galax- 
ies is typically less than the prediction of primordial nucle- 
osynthesis (e.g., Ettori & Fabian 1999; Ettori 2003; He et 
al. 2005; McGaugh 2007). This discrepancy becomes more 
difficult to reconcile in the cores of galaxy clusters (e.g., 
Sand et al. 2003). The best fit of cosmological parame- 
ters with tiny temperature fluctuations of the cosmic mi- 
crowave background (CMB) radiation and large-scale struc- 
ture clustering shows that relative to the critical mass den- 
sity p c in the universe, the mass density of baryon mat- 
ter is n b = 0.0224 ± 0.0009ft M „ and the total matter den- 
sity is fi m — 0.135l ( Q;Qog/i^ 2 , where parameter hum is re- 
lated to the Hubble constant Ho by Ho = 100/tioo km s _1 
Mpc -1 . Therefore, the mean cosmic baryon mass fraction is 
fb = Q, b /Q, m = 0.166±o;ol3 ( e -g-, He et al. 2005 and refer- 
ences therein). While there are different methods in deter- 
mining the ft value, the cosmic baryon fraction fb is around 
0.17 (e.g., McGaugh 2007). However in clusters of galaxies, 
the average gas (baryon) fraction inferred by two methods 
are about 0.107±^ and 0.111±g;g|l (e.g., Ettori 2003). 
Others estimated that the baryon fraction fb observed in 
clusters of galaxies can be lower than the cosmic baryon 
fraction by about 10% - 20% at z = (e.g., He et al. 2005). 
Some even claimed that the value of fb can be lowered by as 
much as 30% (e.g., Ettori 2003). In conclusion, the baryon 
fraction in most clusters of galaxies are systematically lower 
than the average cosmic value fb except those highest es- 
timates for gas (baryon) mass fraction in some clusters of 
galaxies (e.g., A426, A2142, RXJ1350; see Ettori 2003). To 
resolve this important issue, the notion of Warm-Hot Inter- 
galactic Medium (WHIM) has been introduced (e.g., Cen 
& Ostriker 1999, 2006; Ettori 2003). In their opinion, the 
WHIM may actually exist within clusters of galaxies to ac- 
count for the mass of 'missing baryons', yet the WHIM can- 
not be detected at present because it does not emit X-rays. 
These results show that a significant fraction (~ 40% — 50%) 
of the baryon component might be found in the form of 
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WHIM in the temperature range of T ~ 10 5 7 K (e.g., Cen 
& Ostriker 2006). 

As will be discussed in more details, this problem of 
'missing baryons' in our model scenario is a natural conse- 
quence of galaxy cluster winds, be it sustained or sporadic 
or be it stationary or dynamic during the evolution of galaxy 
clusters. These so-called 'missing baryons' are blown away 
in the form of hot ICM and cool down gradually with time; 
with relatively low temperatures, they should mostly reside 
in the periphery of galaxy clusters and spread out in space 
on much larger scales. Meanwhile, the dark matter halo may 
contract within clusters of galaxies in our model. Therefore 
the mass fraction of baryons fb (i.e., the mass ratio of to- 
tal baryons to the total gravitational mass inferred) would 
be lower than the initial value when a cluster of galaxies 
was born and started to evolve. The age of galaxy clus- 
ters is estimated to fall in the range of ~ 10 9 — 10 10 yr 
(e.g., Fabian 1994). As galaxy cluster winds may have ex- 
isted since galaxy clusters were born, the timescale of galaxy 
cluster winds would be comparable to or somewhat less than 
this estimate. In Section 3, we show a few specific examples 
of numerical shock flow solutions in our model and estimate 
the loss of baryons within a timescale of ~ 10 9 yr. 

As different behaviours of temperature profiles have 
been inferred from X-ray observations of galaxy clusters 
(Markevitch 1996 and Markevitch et al. 2005 for decreas- 
ing temperatures with increasing radius; Peres et al. 1998 
and Sanders & Fabian 2006 for nearly constant tempera- 
tures in several galaxy clusters; McNamara et al. 2005 and 
Blanton et al. 2001 for increasing temperatures with increas- 
ing radius) and electron number densities are observed to fit 
a power law fairly well (e.g., Peres et al. 1998; Nulsen et al. 
2005b), we shall take the specific entropy conservation along 
streamlines as the equation of state and see how well this 
may account for the various observed profiles of thermody- 
namic variables. By properly choosing model parameters in 
various regimes, we can describe properties of galaxy clus- 
ters to a considerable extent. 

This paper is structured as follows. The background 
and motivation of our model development is introduced in 
Section 1. Section 2 presents in order the basic formulation 
for the two-fluid model of spherical symmetry involving two 
polytropic fluids, self-similar transformation, asymptotic so- 
lutions at small and large x, singular surfaces and sonic crit- 
ical curves, and shock conditions. In Section 3, we show nu- 
merical examples of quasi-static solutions for three different 
situations. The major results are summarized in Section 4. 
Finally we discuss our model results and numerical solutions 
in Section 5. Certain mathematical details are contained in 
Appendices A through G for the convenience of reference. 



2 MODEL FORMULATION 

As theoretical idealization and simplification, a cluster of 
galaxies is approximated as fully relaxed or virialized and 
usually modelled as a static equilibrium system in radial 
force balance with spherical symmetry. However, large-scale 
shock features as observed in clusters of galaxies reveal the 
presence of large-scale flows, although these implied large- 
scale flows may not be directly measurable at this stage. 
Our main motivation of this model analysis is to provide a 



class of dynamic (rather than static or stationary) models 
for clusters of galaxies with spherical symmetry. We hope to 
understand a few basic aspects of this dynamic model frame- 
work. We have several plausible processes in mind. First, the 
formation of clusters of galaxies through large-scale grav- 
itational collapse involving dark matter and baryon mat- 
ter. Secondly, activities of a central AGN (involving accre- 
tions of baryon matter as well as dark matter; e.g., Hu et 
al. 2004) onto supermassive black holes may give rise to a 
quasi-spherically symmetric component of disturbances on 
large scales which can evolve into shocks. Thirdly, merging 
processes may reach a later phase of core confinement such 
that a large-scale quasi-spherical symmetry may be a sensi- 
ble approximation; while releasing energy, it takes time for 
such a dynamic system to relax and adjust itself. 



2.1 Self-Similar Equations for a Two-Fluid Model 

To study dynamic behaviours of visible baryon matter (such 
as X-ray and radio emissions from the hot ICM) under the 
joint gravity of both massive dark matter and baryon mat- 
ter together, we adopt three assumptions for the dark mat- 
ter halo. First, self-interacting dark matter particle models 
have been proposed earlier by some researchers (e.g., Carlson 
et al. 1992; Machacek 1994; Spergel & Steinhardt 2000) to 
solve the problems encountered by cold dark matter models. 
Furthermore, properties of collisional dark matter particles 
(e.g., Ostriker 2000; Hu & Lou 2007) and fluid dark matter 
(e.g., Peebles 2000; Subramanian 2000; Moore et al. 2000; 
Hennawi & Ostriker 2002; Lou 2005; Hu et al. 2006) were 
proposed as an alternative approach to probe DM dynam- 
ics. On large scales, one may view high velocity dispersions 
(~ 700 — 1000 km s _1 ) of DM particles to produce an effec- 
tive pressure against gravity as described by the Jeans equa- 
tion (e.g., Binney & Tremaine 1987). In particular, Evrard 
et al. (1994) numerically simulated formation of galaxies us- 
ing a model consisting of two gravitationally coupled fluids 
representing dark matter and baryon matter. While a distri- 
bution function approach can be applied to study properties 
of galaxy clusters, we model a dark matter halo in a 'fluid' 
approximation to simplify the mathematical treatment. On 
large scales and without resonances, we should be able to 
understand various dynamic behaviours of hot intracluster 
medium (ICM) and the dark matter halo in this two- fluid ap- 
proximation (Lou 2005). Secondly, we assume the two- fluid 
system of galaxy clusters to be spherically symmetric with a 
common centre and without rotation for simplicity. On much 
smaller scales, this simplification would be a very good ap- 
proximation for globular clusters containing millions of stars 
and gas. Thirdly, dark matter interacts with hot ICM only 
through gravity. Based on the above assumptions, we read- 
ily write out a set of coupled nonlinear partial differential 
equations to describe the two- fluid flow system with hot ICM 
and dark matter halo coupled by gravity. In spherical polar 
coordinates (r, 0, <f>), the equation for mass conservation is 
described by 



dMi dMi 
dt 1 dr 







and 



2 

— = 4.r Pi 



(1) 



where r is radius and t is time; subscripts i = 1, 2 stand 
for dark matter halo (fluid 1) and hot ICM (fluid 2), re- 
spectively. For simplicity, all variables with a subscript i 
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denote associations with fluid i. When referring to a specific 
fluid i, we would write out subscripts 1 or 2 explicitly. Here, 
Mi(r,i), pi(r,t), and m(r,t) are respectively the enclosed 
mass, the mass density, and the radial flow speed at radius 
r and time t for fluid i. Another familiar form equivalent to 
the above continuity equation fTJ is 



dpi 19,2 



: 



(2) 



dui du, 
dt 1 dr 



(3) 



The Euler radial momentum equation is 

1 dPj _ G(Mi + M 2 ) 
pi dr r 2 

where Pi(r, t) is the pressure at radius r and time t for fluid 
i; Mi and M 2 are the enclosed masses of fluid 1 and fluid 
2, respectively; G = 6.67 x 10~ n kg _1 m 3 s~ 2 is the gravi- 
tational constant Q The coupling effect of the two fluids is 
explicitly contained in the gravity term in equation ([3]). Fi- 
nally, we take the conservation equation of specific 'entropy' 
along streamlines as the equation of state (e.g., Fatuzzo et 
al. 2004; Wang & Lou 2007; Lou & Cao 2007) 



d d\ ( Pi . 



(4) 



where 7; is the polytropic index for fluid i. In general, poly- 
tropic indices 71 and 72 are allowed to be different. 

Equations (HJ — (|4j) form a set of nonlinear partial dif- 
ferential equations and contain an important subset of non- 
linear self-similar solutions with or without shocks. We now 
introduce a set of self-similar transformation for the two 
polytropic fluids (e.g., Suto & Silk 1988; Lou & Wang 2006) 
below 



Pi 



ai(xj) 
A-KGt 2 



Mi = — trti(ii) 



(5) 



(3n» - 2)G 

where x% is the independent dimensionless similarity variable 
for fluid i (x\ and X2 are actually related to each other); Ki 
is a scale parameter of self-similar transformation for fluid i; 
Hi is another index parameter of self-similar transformation 
for fluid i noted above; ou(xi) is the reduced mass density 
for fluid i; Vi(xi) is the reduced radial flow speed for fluid i; 
Pi(xi) is the reduced pressure for fluid i; and rm(xi) is the 
reduced enclosed mass for fluid i. All these reduced variables 
are functions of independent variable Xi only. In principle, 
parameter m does not need to be equal n 2 . But equation 
([3]) contains factors like t ni and t n2 at the same time after 
the self-similar transformation. Therefore in order to obtain 
the self-similar dimensionless ordinary differential equations 
(ODEs) in terms of Xi without involving explicit temporal 
factors of time t, we simply require n\ — 722 = n as a single 
parameter and shall not distinguish the two from now on. 
With the specific entropy conservation Q along streamlines, 
it is not necessary to require n = 2— 74 here (e.g., Yahil 1983; 
Suto & Silk 1988; Lou & Wang 2006, 2007). In fact with 1 < 
7i < 2, we shall explore three possible situations of 2/3 < 
n < l,n= 1 and n > 1, respectively. We note that self- 
similar processes of galaxy cluster evolution and of cooling 

1 In Lou & Wang (2006), there is a typo in the unit of G; it 
should be cm 3 instead of cm -3 . 



waves in galaxy clusters have also been studied previously 
(e.g., Bertschinger 1989; Jain & Bertschinger 1996). Here, 
cooling waves refer to a self-similar expansion of cooling flow 
region. While the cooling region expands, the hot ICM itself 
does not move out. 

By performing self-similar transformation the re- 
duced enclosed mass rrii(xi) can be expressed as (e.g., Suto 
& Silk 1988; Lou & Wang 2006; Lou & Cao 2007) 



rrii{xi) = ctiXi (nxi - Vi) 



(6) 



Since the enclosed mass Mi > 0, we should require > 
and nxi > Vi for n > 2/3, while for n < 2/3, we require 
rrii < and nxi < Vi\ the latter is generally impossible for 
semi-complete solutions in the range of + < x < +00. Us- 
ing equation ([SJl and self-similar transformation (JSJ), specific 
entropy conservation along streamlines leads to 



Co,imf(xi) , 



(7) 



PijXi) 

af (xi) 

where qi = 2(n + 7* — 2)/(3n — 2) is a naturally emerged 
index parameter and Co,i is an integration constant for each 
fluid i. For 7, 7^ 4/3, we can always effectively combine the 
two coefficients Co,j and Ki into a new single constant coef- 
ficient, corresponding to a coefficient rescaling in self-similar 
transformation ([5]). The case of a single fluid with 7 = 4/3 
is separately considered by Lou & Cao (2007; see also Gol- 
dreich & Weber 1980 and Yahil 1983). It then suffices to 
consider the equation of state in the form of 

Pi{Xi) 



&i (Xi) 



= mf{xi) 



for fluid i with coefficients Co,i being absorbed without loss 
of generality. Apparently, there are two linearly related inde- 
pendent variables x\ and X2 respectively for the two coupled 
fluids under consideration. We introduce a convenient ratio 
k = (-ft'i/-K 2 ) 1//2 such that k = xijx\. We can then express 
all dependent variables as functions of x\ only. From now 
on, we shall rewrite xi as x for simplicity and thus X2 = kx. 
Now mass and momentum conservation equations ((2]) and 
((3| can be straightforwardly cast into the following dimen- 
sionless forms of four ODEs, namely 



, .dai dvi n (x-vi) 

(nx - Vi)— ai—— = -2 

ax ax x 

(n — l)vi — (nx — Vi)— 7— = 

ax 



Q-l 



(9) 



■ 7i« 



<Jl+71-2„2<7i 



1 \gi dati 

(nx — vi) — — 
ax 



giaf +11 ' 1 x 2qi (nx - Ui) ?1_1 (3n - 2) 
ai(nx — vi) a 2 (Knx — « 2 ) 



(3n - 2) 



(3n - 2)k 



, , dd2 d,V2 „(kx — V2) 

(nnx — v 2 )— ol%—— = —2 a 2 , 

ax ax x 



(10) 
(11) 



(n — l)v2 



(nnx — V2) dv2 
k dx 



72 a 92+72-2 K 2g 2a ,2g 2| 



nnx — V2 ) 



, dot2 
dx 



q2al 2+12 - 1 K q2 x 2q2 (nnx - V2) q2 ~ 1 (3n - 2) 
a2(nnx — v 2 ) kol\(tix — v±) 



(3n - 2) 



(3n - 2) 



(12) 
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After rearrangement and algebraic manipulations of above 
four equations (|9} — (|12[) > the reduced radial flow speeds v\ (x) 
and V2(x) and the reduced mass densities ai(x) and 02(2:) 
are then determined equivalently by four coupled first-order 
nonlinear ODEs shown below. 

dm (as) _ .Ai(a;) 



(14) 
(15) 
(16) 



dx 




dvi(x) 


Vi(x) 


dx 


T>i{x) 


dct2(x) 


A 2 (x) 


dx 


T> 2 (x) 


dV2(x) 


V 2 (x) 


dx 


V 2 {x) 



where explicit expressions of denominators D\ (x) and T>2 (x) 
and numerators Vi(x) and V2(x) are shown in Appendix 
A. This set of nonlinear ODEs contains various types of 
asymptotic solutions at large and small x. We mainly focus 
on two of them, one at large x and the other at small x. 

2.2 Tests of Our Model Formulation 

Our model formulation here is sufficiently general and can 
be readily reduced to various known formulations through 
various paths of reduction. We indicate these different paths 
of reduction below to test and confirm the robustness of our 
approach. 

For n = 1 and 71 = 72 = 1, the formulation here re- 
duces to a spherical composite system of two isothermal flu- 
ids coupled by gravity as explored by Lou (2005) . Static and 
dynamic models may be constructed in that relatively sim- 
ple theoretical framework in contexts of clusters of galaxies 
or of globular clusters but on much smaller scales. 

By setting physical variables of one of the two fluids 
to vanish, the formulation here reduces to that of a single 
spherical system for a more general polytropic fluid under 
self-gravity (e.g., Fatuzzo et al. 2004; Lou & Cao 2007; Wang 
& Lou 2007 in preparation). 

If we further set n+7 = 2 for a single general polytropic 
gas under self-gravity, our formulation here reduces to that 
of a conventional polytropic gas with a constant specific en- 
tropy distribution in time and space (Goldreich & Weber 
1980; Yahil 1984; Suto & Silk 1988; Lou & Gao 2006; Lou 
& Wang 2006, 2007). 

The single isothermal fluid of spherical symmetry corre- 
sponds to n — 1 and 7 = 1 in our formulation. This problem 
has been extensively explored in the literature (e.g., Larson 
1969a, b; Penston 1969a, b; Shu 1977; Whitworth & Sum- 
mers 1985; Hunter 1977, 1986; Tsai & Hsu 1995; Shu et al. 
2002; Lou & Shen 2004; Shen & Lou 2004; Fatuzzo et al. 
2004; Bian & Lou 2005; Yu & Lou 2005; Yu et al. 2006). 

Our model is formulated for two gravitationally cou- 
pled fluids which are polytropic in the more general sense 
(see equation |3J); this holds the key difference in reference to 
previous single fluid model under self-gravity with various 
equations of state. We will show in this section by examples 
that the solutions similar to the well-known solutions in the 
single fluid framework (such as the static solution, the cen- 
tral free-fall solution (Shu 1977) and the Larson— Penston 



type solution (Larson 1969a, b; Penston 1969a, b)) can also 
be derived within our model framework. 

The static solution in which both flow velocities of the 
hot gas and dark matter vanish throughout the entire space 
can be found in our model framework. This is an exact 
global solution similar to the static solutions with central 
divergence (i.e., singular isothermal sphere (SIS) and singu- 
lar polytropic sphere (SPS)) in single fluid model framework 
(e.g., Shu 1977; Cheng 1978; Lou & Shen 2004; Lou & Wang 
2006; Lou & Cao 2007). In our two-fluid model framework, 
the global static solution (i.e., singular double polytropic 
spheres (SDPS)) is simply 



Vi = V2 = , Qi = A\X 



-2/r, 



OL2 — A2X 



-2/r) 



(17) 



Here, the two positive density coefficients Ai and A2 for two 
static fluids are readily determined by the following pair of 
equations 



n{M + A 2 ) 
2(3n - 2) 



{2-n)n qi - L A\ 



91-1 Aqi+H- 



njAi + A 2 ) _ 3 92 -2, 9 
2(3n-2) 1 



92-1^92+72-1 



(18) 



where qi = 2(n + 74 — 2)/(3n — 2) and the ratio k for the 
two fluids breaks the symmetry of the above two relations; 
this symmetry would be explicit for n — 1. For physical so- 
lutions with positive Ai and A2, it is necessary to require 
2/3 < n < 2. In other words, both power-law density scal- 
ings fall between x~ x and x~ 3 . For a set of four specified 
parameters (n, 71, 72, k), equation (|18p does possess sen- 
sible real solutions for both coefficients Ai > and A2 > 
and our model then gives a singular static solution for both 
polytropic fluids simultaneously, i.e., two gravity coupled 
singular polytropic spheres (SPSs) with divergence of mass 
densities as x — » + . In astrophysical applications, we need 
to introduce a proper central cutoff. 

The LP asymptotic solution and Shu's central free-fall 
asymptotic solution were constructed for a single isothermal 
sphere; in order to get analogous asymptotic solutions in 
our two-fluid model, we set scaling index parameter n = 1 
(defined in self-similar transformation equation (0) and the 
polytropic indices of both fluids 71 = 1 and 72 = 1. This is 
a self-gravitating system of two coupled singular isothermal 
spheres (Lou 2005). 

For two coupled isothermal gas spheres with « = 1, the 
generalized version of the Einstein-de Sitter solution (Whit- 
worth & Summers 1985; Shu et al. 2002; Lou & Shen 2004; 
Lou & Zhai 2007 in preparation) in our model is 



vi — V2 — 2x/3 , Qi = a.2 — 1/3 



(19) 



This isothermal solution is an exact global solution. For k 7^ 
1, this kind of solutions does not exist. 

For two coupled conventional polytropic spheres with 
n + 7i=n + 72 = 2 and « = 1, the generalized version of 
the Einstein-de Sitter solution of our model is still given by 
equation JT9J (Lou & Wang 2006; Wang & Lou 2007; Lou 
& Cao 2007). 

The asymptotic behaviour when x — > + of the solution 
similar to Shu's central free-fall asymptotic solution in our 
model framework is 



Ttix 



-1/2 



v 2 



H2X 



-1/2 



(20) 
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on — » Qix 



-3/2 



Ol2 — » 62a; 



-3/2 



(21) 



Here, the two relations among the four coefficients TLi, H.2, 
Qi > and C/2 > are 



H? + 2fti£i + =0 , 



+ 2«W 2 & + I^TiiQi = 



(22) 



(23) 



Therefore, there are only two free coefficients for this kind of 
central free-fall asymptotic solutions as x — > + . The mass 
at the centre for fluid 1 are determined by mi(0) = —TiiQi 
and the mass at the centre for fluid 2 are determined by 
ti 2 (0) = —k^HiQi- Therefore, Hi and TI2 are both negative 
in equation (|20p . corresponding to central free-falls of both 
fluids. 

For a polytropic single fluid model, quasi-static asymp- 
totic solutions have been constructed by Lou & Wang 
(2006). In the next subsection, we will also construct similar 
quasi-static asymptotic solutions in our two-fluid model and 
apply the model to the large-scale dynamics of hot gas and 
dark matter in clusters of galaxies in this paper. 



2.3 Construction of Two-Fluid Quasi-Static 
Solutions 

As the sound speed of hot ICM in clusters of galaxies is 
of the order of ~ 1000 km s _1 and may decrease towards 
the centre (such as Abell 2052 in Blanton et al. 2001, and 
MS0735. 6+7421 in McNamara et al. 2005 and in Cavaliere & 
Fusco-Femiano 1978) while the observed radial flow speed of 
hot ICM close to the centre is not large, we then construct 
the quasi-static solution from equations fll3f) — (|16[) in the 
regime of small x (i.e., parallel to the quasi-static solution 
of a single fluid of Lou & Wang 2006 and Wang & Lou 
2007)0 

Now we take the static SPS solution (|17|l and (|18[) as 
given in the subsection 12.21 as the leading-order term of the 
asymptotic quasi-static solution as x — ► + . Expanding to 
the second order of this series, the quasi-static solution for 
two gravity coupled polytropic fluids then takes the form of 



Vi — * L\x + . 
v 2 — > L 2 x k + . 



(24) 



A lX ~ 2/n + N lX R + 
^' 2/n + N 2 x R + 



Ql 

a 2 A 2 x' 2/n + N 2 x R + ... , (25) 

where R = k — 1 — 2/n and index parametei0 k is determined 
by the following quartic algebraic equation 

(26) 



C M fc 4 + C fe ,2 A; 3 + C fc ,3 k 2 + C kA k + C k , 5 = 



Explicit expressions of the five coefficients {Ck,i, Ck,2, Ck,3, 
Cfe,4, Ck,s) can be found in Appendix B. In general, the 
roots of algebraic quartic equation (|26[1 can be complex. For 
the higher order terms with respect to static SPS solution, 
we should require Re(k) > 1 (in static SPS solution (I17f) . 

2 Yahil & Ostriker (1973) discussed a steady outflow of gas from 
galaxy clusters. Their steady wind results show that the gas ve- 
locity towards the centre is also very small (see their Fig. 3). 

3 The index parameter k is taken to be the same for both fluids. 



vi = v 2 = Ox is regarded as the zeroth order); and the two 
coefficients L\ and L 2 are two free parameters with the other 
two coefficients Ni and N 2 readily determined by 



Ni = 



N 2 



[2(n - l)+nk]A 1 L 1 

(k - l)n 2 
[2(n- 1) +nk]A 2 L2 

n(k — l)n 2 



(27) 



For a complex index parameter k, please refer to Appendix 
C for quasi-static solutions with asymptotic oscillations in 
the regime of small x (see also Lou & Wang 2006 for quasi- 
static solutions with asymptotic oscillatory behaviours). 



2.4 Asymptotic Solutions at Large x 

In this section, we shall derive asymptotic behaviours of ra- 
dial flow speed and mass density of the two fluids from equa- 
tions (|13|) — (|16|) as x — » +00. 

The asymptotic series solution at large x takes the form 

of 



ai(x) -> E lX - 2/n + hx~ 3/n + . 
a 2 (x) -> E 2 x~ 2/n + I 2 x- 3/n + . 



v\(x) — » H\x 



V2(x) — ► H2X 



-1/n+l 



-1/n+l 



+ GlX 



+ G 2 X 



-2/n+l 



-2/n+l 



+ . 



+ . 



(28) 
(29) 
(30) 
(31) 



where 2/3 < n < 2 is required and the four coefficients 
Ei, E2, Hi and H2 are fairly arbitrary, while the other four 
coefficients Ii, I2, Gi and G2 can be expressed in terms of 
these four arbitrary coefficients (details of these expressions 
are contained in Appendix D). 

For asymptotic radial flow speed solutions (|30[) and 
(|31|) . the flow velocities diverge asi-> +00 for n > 1 unless 
Hi = H2 = 0. For a real astrophysical system, its size is 
finite; and we may need to introduce a spatial cutoff at a 
given time in order to make use of these solutions for n > 1 
with Hi 7^ and H2 7^ 0. For example, the typical size of a 
galaxy cluster is of the order of several to ten Mpcs. 

We are certainly interested in finite asymptotic solu- 
tions J30j and ||3TJ with n = 1 and 2/3 < n < 1. While for 
n > 1, it is possible to set Hi = H2 = 0, and the asymptotic 
solutions finite at large x become 

Eix-' 2/n + Fix- 4/n+1 



ai(x) 



a 2 {x) -> E 2 x~ 2/n + F 2 x 



-4/n 



+ ■ 



+ ■ 



vi(x) — > Gix 



V2(x) — > G2X 



-2/n+l 



+ Dix 



-4/n+2 



+ 



-2/r, 



_1 + D 2 X 



-4/n+2 



(32) 



(33) 



(34) 



(35) 



where Ei and E2 are two fairly arbitrary constants and the 
condition n < 2 is still required; Gi and G2 are still defined 
in Appendix D. Fi, F2, Di, and D2 are four constant coef- 
ficients of next order expansion terms; these coefficients are 
determined by specified values of Ei and E2 (further details 
of these coefficient expressions can be found in Appendix 
D). In principle, we can carry out this series expansion for 
large x to the desired order if needed. 
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2.5 Singular Surfaces and 
Sonic Critical Curves 

The singular surface and sonic critical curve in a single gas 
flow have been analyzed in details by previous authors (e.g., 
Suto & Silk 1988; Whitworth & Summers 1985; Lou & Wang 
2006; Wang & Lou 2007). For a single gas flow of spherical 
symmetry, there are smooth solutions going across the sonic 
critical curve analytically (e.g., Suto & Silk 1988; Whitworth 
& Summers 1985; Lou & Shen 2004) or with shocks (e.g., 
Tsai & Hsu 1995; Shu et al. 2000; Bian & Lou 2005; Lou & 
Wang 2006; Yu et al. 2006; Wang & Lou 2007; Lou & Cao 
2007). In our two- fluid model, each fluid component has its 
own singular surface and sonic critical curve; they are not 
the same for the two fluids in general and should be treated 
separately. 

The singular surfaces are a set of points when T>i(x) = 
for fluid 1 by definition (fXo) and V 2 (x) = for fluid 
2 by definition (|A6p . Physically, when the travel speed of 
disturbances relative to the flow speed is equal to the local 
sound speed in a fluid, we encounter a singular surface. As 
there are two fluids in our model, there are two singular 
surfaces possible for a given set of relevant parameters. The 
singular surface {x, Vi, ai} of fluid 1 as defined by T>i(x) = 
is then given the following equation 



(nx-v-i) --yia\ 



91+71- 1 ^, 2 9i 



x 91 (nx — vi) q 







(36) 



and the singular surface {x, V2, a 2 } of fluid 2 as defined by 
T^2(x) — is given by a similar equation 

(kux — V2)' 2 — K 292 72a| 2+72_1 a; 2 ' 32 (nnx — U2) 92 = . (37) 

Mathematically, sonic critical curves are characteristic 
lines on the singular surface with the numerators and de- 
nominators of self-similar nonlinear ODEs (|13|l and (I14p . 
and (|15l) and (|16p separately being zero simultaneously. The 
intersection of such two surfaces in {x ,Vi ,a{\ space gives 
the sonic critical curve in each fluid i. As noted above, when 
the travel speed of disturbances relative to the local flow 
speed of a specific fluid is equal to the local sound speed, a 
singularity arises. If one wants to get continuous flow solu- 
tions throughout the entire range of x, such solutions must 
pass across the singular surface via the sonic critical curves. 

In our two-fluid model, as Vi(x) and Ai(x) contain a 2 
and v 2 , we need proper values of 012 and V2 in order to de- 
termine the sonic critical curve of fluid 1; in reciprocal and 
in parallel, as V2(x) and A2(x) contain ai and Vi, we need 
proper values of a\ and v\ in order to determine the sonic 
critical curve of fluid 2. This is a major yet expected differ- 
ence as compared with the case of a single fluid. Therefore 
in principle, for a given position of x, there is a critical point 
(a point on the critical curve) for fluid i once a set of (an, 
Vi) at that x of the other fluid is given. As the two fluids do 
not have their critical curves at the same x in general, there 
are only two eigendirections for a given set of parameters, 
which is similar to the case of a single fluid model (Lou & 
Wang 2006; Wang & Lou 2007). The equations and specific 
procedure to determine the eigendirections can be found in 
Appendix E. 

In reference to our quasi-static solution for small x, the 
condition for the corresponding singular surface T>\(x) = 



of fluid 1 has the asymptotic behaviour when x — + + : 

/ !l+7l -i 2 <?1 \ 1/(2 ~ 91) 
nx = [7iaj x I +vi 

_ J 7iA?1 +-, 1 -lj 1/(2 -' l) ;c l-2/[n(2- gi )] _ 



(38) 



For T>2(x) — 0, the asymptotic behaviour of the corre- 
sponding singular surface for fluid 2 when x — > + is 

1/(2-92) 



K/,'.r = I K 2 ' 72 72a| 2+72 



+ V2 



292 492+72-11 1/(2 92 ' l-2/[n(2-9 2 )] 



72^ 2 ^ 



(39) 



Therefore if 1 — 2/[n(2 — qt)] > 0, the corresponding singular 
surface of fluid i passes through the origin point (x = and 
Vi =0); otherwise, it cannot pass the zero point. 

When x —* +00, the asymptotic behaviour of cti(x) is 
(see asymptotic solutions (J25J, ([29]), (|D9|| and (|D10]|) char- 
acterized by ai — > EiX~ 2 ^" + . . .. Then for this same limit of 
large x, the singular surface condition T>\(x) = becomes 



vi = nx — 



l/(2-9l) 



- nx-^El^-^^x 1 - 2 ^ 2 -^ . (40) 

Since when x — > +00, v\(x) is either zero or Vi — > Hix 1 ^ 1 ^ 71 
and 1 — 2/[n(2 — qi)] 7^ 1 for any values of n and qi, our 
asymptotic solutions will not lie on the singular surface for 
fluid 1. For T>2(x) — 0, we have the similar result 

v 2 = nnx - ( K 292 72a!f +72 -V 92 ) 1/<2 ~ 92) 

^ K na;- (72 K 292 J B2 92+72 ~ 1 ) 1/(2 ~ 92 ^ 1 - 2/[n(2 - 92)1 . (41) 

Since when x — » +00, 112 is either zero or V2 — » H2X X - 1/n and 
1 — 2/[n(2 — 32)] 7^ 1 for any values of n and 52, our asymp- 
totic solutions will not lie on the singular surface for fluid 
2. In summary, for either fluid in the regime of x — * +00, 
our asymptotic solutions will not encounter the respective 
singular surfaces. 



2.6 Jump Conditions for Self-Similar Shocks 

Large-scale shocks have been revealed in clusters of galax- 
ies through high-resolution X-ray imaging observations and 
radio observations (e.g., Nulsen et al. 2005a, 2005b; McNa- 
mara et al. 2005; Bagchi et al. 2006). In order to probe and 
model such shock features in clusters of galaxies, we con- 
struct numerically semi-complete flow solutions with shocks 
across singular surfaces in both hot gas and dark matter 
halo separately. These large-scale shocks travel outward in a 
self-similar manner with asymptotic flow signatures at large 
r for given time t. 

In principle, the self-similar form of shock radial posi- 
tions r 3; i in fluid i can be expressed either by the down- 
stream parameters (Ki and Xi on the downstream side of 
a shock) or by the upstream parameters (Ki and Xi on the 
upstream side of a shock). For simplicity, we express shock 
radial positions r Si ; in terms of the upstream parameters. 
For fluid i, the shock radial position is denoted by 

'1/2, 



- K ''t'W 



(42) 



where K U: i is the upstream value of Ki and x Ut i is the shock 
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position in terms of the upstream independent self-similar 
variable for fluid i. We introduce the sound speed ratio Ai 
across a shock for fluid i as 



Ud,i = \i~K]Ht n 1 v d> i(x d , l ) , 



A, = 



Kd 



1/2 



(43) 



where the subscript d denotes the downstream side (i.e., 
from the shock position towards the centre) while the sub- 
script u denotes the upstream side (i.e., from the shock po- 
sition towards infinity). In reference to a shock, Kd,i is the 
parameter of Ki on the downstream side; Xd,i is the shock 
location in terms of the downstream self-similar variable for 
fluid i; pd,i is the mass density downstream of a shock in 
fluid i while a d ,i is the reduced mass density of a shock 
downstream of fluid i; Pd ; i is the pressure downstream of 
a shock for fluid i while f3 d ,i is the reduced pressure down- 
stream of a shock for fluid i; Ud,i is the radial flow speed 
downstream of a shock for fluid i while Vd,% is the reduced 
radial flow speed downstream of a shock for fluid i; Md,i is 
the enclosed mass downstream of a shock for fluid i while 
m d>i is the reduced enclosed mass downstream of a shock for 
fluid i. In one-to-one correspondence, all these variables with 
subscript {u, i} refer to variables on the upstream side of a 
shock for fluid i. In order to construct self-similar shocks, we 
require index parameter n to be the same across a shock to 
avoid unphysical interface separation. As the shock radius is 
r s ,i = K Uti t n x Ui i, the dimensional shock speed is then given 
by 



dr s 



»0" 



m .--„,. -.« = n-Ji. (44) 

Here, u s ,i is the travel speed of a self-similar shock in fluid 
i, indicating that the shock actually travels with a variable 
speed for n ^ 1. For n > 1 and n < 1, the shock travel 
speed increases (acceleration) and decreases (deceleration) 
with time t, respectively. For n = 1 and fairly arbitrary ji, 
such a self-similar shock (not necessarily isothermal though) 
travels with a constant speed. In clusters of galaxies, posi- 
tions of shocks observed vary from tens of kpcs (e.g., McNa- 
mara et al. 2005; Nulsen et al. 2005b) to a few Mpcs (e.g., 
A3667 in Rotteringer et al. 1997; galaxy cluster Abell 3376 in 
Bagchi et al. 2006; ZwCl 2341.1+0000 in Bagchi et al. 2002; 
A786, A2255, A2256 in Ensslin et al. 1998). In our scenario, 
such shocks actually travel to current positions from a in- 
ner region around the cluster centre after their emergence, 
which can be estimated in our model framework. Then in 
our scenario at time t, the shock position r s>i is determined 
by the following equation 



(45) 



1 /2 

where C s % = K ' ■ X u ,i IS cl constant to be estimated from 
observations at a certain time under specific situations. 

The self-similar transformation for variables on the 
downstream side of a shock is 



Xd,i — 



Pd,i = 



a d ,i(x dii ) 



Pd,i = Aj JV " 4 'J G Pd,i{ x d,i) 



(46) 

(47) 
(48) 



M d ,; 



A? 



K 3 /h 3n ~ 



-m d ,i{x d ,i) 



(49) 



(50) 



(3n - 2)G 

In the shock reference framework, the mass conservation 
equation across a shock front is 

pd,i(u s ,i - u d ^) - p u ,i(u s ,i - u u ,i) = , (51) 
which can be rewritten conveniently as 
[Pi(u s ,i - Ui )]t = ; (52) 
likewise, the radial momentum conservation then gives 



[Pi + Pi(u s ,i 







and the energy conservation equation leads to 



Pi(u s ,i - Ui) 



= 



(53) 



(54) 



where the pair of brackets denotes the difference of the argu- 
ment on the downstream (supercript d) and upstream (sub- 
script u) sides of a shock. As for a single fluid, we now in- 
troduce two new variables Yd,i and F Uii below 

r d ,i = n - v d ,i/x d ,i , (55) 

T u ,i =n — v u ,i/x u ,i ■ (56) 

Once we know the values of (F, a, x) on the downstream side 
(indicated by a subscript d) of a shock, we can immediately 
calculate the corresponding variables on the upstream side 
(indicated by a subscript u) or vice versa. Details of shock 
calculations can be found in Appendix F. We only show the 
major results here. The variable F Ui i on the upstream side 
can be computed from the variables on the downstream side 
from the following equation 



2 7l 



,,+ 7 ,-l r ,,-l 3,,-2 {li - 1) p 
. u d,i 1 d,i x d,i T i i \ , 



( 7i +i)-<w 1 ( 7i +ir°" ' (57) 

other variables can be readily determined in a straightfor- 
ward manner 

&d,i^d ,i /ro\ 

a u ,i = — , (58) 



<7,+7» r <?, r 3<7,-2 , _p2 

- 1 - A a^ri i I ^d^L d,l 



d,i 
Old 



,i^u,i^d,i^ 



d^dji 

1/(394-2) 



l/(39i-2) 



It then follows that 

^u,i — Xu^iiri r^ ; j) 



.(59) 



(60) 



and the ratio Ai across a shock in fluid % can be determined 
accordingly. The upstream Mach number M u ,i is defined by 

1/2 

Mu.t = =1 — ^ — I (««,»-««,») 



J- U.l 



1/2 (9i+7i-l)/2 x 1i 



(61) 



where the polytropic sound speed on the upstream 

side of a shock in fluid i, namely 

<>u - \ -i : ( iif^y /2 . (62) 



( dPu,i 

\dp u ,i 



Pu,i 
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In general, shock positions in the two fluids are different 
and the ratio k will change once either fluid goes across the 
singular surface via a shock. Right at the shock front, we 
have r dtl = r U)i = r a ,i, X t = {K d ^/ K U)i ) 1/2 = x u ^/x d>i ; we 
can then recalculate k with the definition k = (K\j K^) 1 ^ 2 = 
X2/X1 where Ki should take the local value. 

The specific 'entropy' of fluid i is conserved along 
streamlines and is given by 

Si = c v ,i In ^~t7^ 

= C, i ln|jrf- 37i)/(3n - 2) (63) 

x(4 7 rG) 7 ^ 1 [(3n-2)G]' Jl Aff I , 

where c Vy i is the specific heat capacity at constant volume for 
fluid i. As the specific entropy increases from upstream side 
to downstream side across a shock, either parameter Ki in- 
creases from upstream side to downstream side for 7^ < 4/3, 
or parameter Ki decreases from upstream side to down- 
stream side for 7; > 4/3; we take ji 7^ 4/3 in this paper. 
For a single self-gravitating polytropic gas with 7 = 4/3, the 
reader is referred to Lou & Cao (2007) for a further theoret- 
ical development of earlier analyses by Goldreich & Weber 
(1980) and Yahil (1983). 

3 RESULTS OF NUMERICAL EXAMPLES 

Up to now, within the self-similar dynamic model framework 
of two gravity coupled polytropic fluids, we have successfully 
constructed the generalized version of quasi-static asymp- 
totic solutions for small x in reference to the model analysis 
of Lou & Wang (2006, 2007). This type of self-similar evo- 
lution eventually approaches a static configuration with a 
diverging density towards to the central core region. Mean- 
while, we have determined the two singular surfaces and the 
shock conditions across the two singular surfaces, respec- 
tively. In order to construct a global semi-complete quasi- 
static solution, the relevant parameters required to be known 
are: the scaling index n as introduced in self-similar transfor- 
mation ©; a proper starting value Xi n i in the small x regime 
to guarantee a reliable numerical integration; the two poly- 
tropic indices of dark matter 71 and of hot ICM 72 respec- 
tively; an estimate of time t when a self-similar evolution is 
presumed to begin; the two parameters K\ and K2 are re- 
lated to the sound speeds of two fluids respectively and are 
introduced in self-similar transformation ([SJ) or equivalently, 
the parameter K2 and the ratio k = (K\f 'A'2) 1 ^ 2 of the two 
fluids; the initial parameters for the quasi-static velocities 
L\ and L2 as defined in the quasi-static solution (|24|l ; the 
independent self-similar variables on the downstream sides 
of shock positions for dark matter x d ,i and for hot ICM x d ,2- 
Once these eleven parameters are specified, a semi-complete 
numerical solution can be established. If we just construct 
a dimensionless solution, then values of K2 and t are not 
needed, indicating that only nine dimensionless parameters 
are required to be known. Note that parameter Xmi needs 
to be carefully chosen. Flow parameters at large x can be 
determined accordingly. 

Now we try to use our two-fluid model to explore dy- 
namic behaviours of hot ICM and dark matter halo for clus- 



ters of galaxies. Throughout this paper, we take fluid 1 to 
represent the dark matter halo and fluid 2 to represent the 
hot fully ionized ICM. Any variables with subscript 2 are 
associated with the hot ICM in galaxy clusters in the model 
analysis. We take the quasi-static solution of both fluids as 
x — » 0+ and go across the singular surfaces with shocks 
in hot gas and in dark matter halo (at different locations 
and thus different outward shock travel speeds) respectively. 
The cluster-scale shocks have been observed in many clus- 
ters of galaxies, which may be related to cluster formation 
processes, central AGN activities or merging of galaxies. Al- 
though shocks in a dark matter halo have not yet been de- 
tected, there is no obvious reason to rule out this possibil- 
ity. Observationally, it may be possible to test their presence 
once the density jump profile of a dark matter halo can be in- 
ferred through effects of gravitational lensing. Very recently, 
Onemli & Sikivie (2007) proposed to interpret certain grav- 
itational lensing observations to be "caustics" (i.e., sharp 
rises of density in DM halos) in galaxy clusters. While dark 
matter shocks in our model differ from such DM caustics 
discussed in the literature, they do share certain similar fea- 
tures and therefore, shocks may also be detected by utilizing 
gravitational lensing effects. More detailed discussion on DM 
caustics can be found in the Discussion section at the end. 
In our model, we can describe various dynamical behaviours 
of hot ICM when r — » +00, including inflow, outflow and 
static solutions. Especially for the outflow solutions at large 
x, we shall refer to them as galaxy cluster winds, just like 
solar and stellar winds or galactic winds on much smaller 
yet different scales. On the basis of galaxy cluster winds 
and flows of dark matter halo, we will estimate the loss of 
baryon matter during a timescale of the order of ~ 10 9 yr 
for the evolution of galaxy clusters. 

As scaling index n is a key parameter controlling asymp- 
totic scaling behaviours of self-similar dynamic solutions, we 
shall discuss model solutions for three cases of 2/3 < n < 1, 
n = 1 and n > 1, respectively, all with 7; in the range 
of 1 < 7, < 2 (polytropic indices 71 and 72 of two fluids 
are allowed to be different in general). We use the standard 
fourth-order Runge-Kutta numerical scheme (e.g., Press et 
al. 1986) to integrate coupled nonlinear ODEs fl"3)l. (fl4|) . 
(| 1 5|) and ()16p with relevant asymptotic solutions in both 
regimes of large and small x. As the analytic quasi-static 
solution is adopted for x — *■ + in this model consideration, 
we start the numerical integration with an assigned value 
of x^i, which is small enoughQ The corresponding initial 
values of on(x) and Vi(x) are determined by quasi-static so- 
lution (|25|) and (|24p with mutually consistent coefficients. 
Once a fluid component encounters its singular surface, we 
let the fluid go across it via a shock; there is a certain de- 
gree of freedom in choosing shock location and thus shock 
speed in constructing solutions. If the solution can be inte- 
grated to infinity (a sufficiently large x in practice), a semi- 
complete solution is then obtained in our two-fluid model. 
Note that in order to avoid unstable numerical integration 
from small x outwards, we need to take some appropriate 
initial parameters and carefully check characteristic features 
of a true quasi-static solution at small x. 



4 The rule of thumb criterion for a sensible choice of Xi n i is that 
Vi(x)/x k remains constant for a certain range of small x. 
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In applications to clusters of galaxies, the virial radius 
corresponds to the radius where the mass density of dark 
matter is 200 times of the critical density p c in the Universe 
(e.g., Navarro et al. 1996). We follow this convention and de- 
note the virial radius by 7-200- Here, the critical mass density 
p c is defined by 



Pc = 3flS/(87rG) 



(64) 



where Ho is the Hubble constant and G = 6.67 x 
HT 11 kg" 1 m 3 s 2 is the gravitational constant. Recent 
Wilkinson Microwave Anisotropy Probe (WMAP) results 
showed a Hubble constant Ho = 72 ±8 km s _1 Mpc -1 (e.g., 
Freedman et al. 2001; Spergel et al. 2003). Thus the critical 
mass density is currently inferred to be p c = 9.7 x 10 -27 kg 
m -3 . In our discussion on clusters of galaxies, the radius r 
will be expressed in unit of the virial radius r2oo- 



3.1 Ranges of K2 and t Values 

Once valueijf] of K2 and time t are chosen in our two-fluid 
model, we can then simulate a dynamic shock flow of hot 
ICM in clusters of galaxies under the gravitational influence 
of a flowing dark matter halo. Physically, we may take t to 
be the timescale when a large-scale shock has emerged. As 
shocks in clusters of galaxies may be initiated by gravita- 
tional core collapses or AGN activities or merging galaxies 
and the timescale for the recurrence of AGNs is of the order 
of ~ 10 8 yr (or equivalently ~2x 10 15 to 2 x 10 16 s, e.g., 
Fabian 1994), we may choose time t of this magnitude order 
in our model applications to clusters of galaxies. Accord- 
ing to dimensionless equation of state (|SJ) and self-similar 
transformation ([5J, we have 

-,2/(392-2) 



K 2 = 



12 

^-M| 2 (47r) 72 - 1 G 92+72 - 1 (3n - 2) 92 



(65) 



For a hot ICM, we simply apply the ideal gas law P2 = 
A/2&BT2, where A/2 is the particle number density, fcs is the 
Boltzmann constant, P2 is the thermal pressure of ICM, and 
T2 is the ICM temperature. Then equation (|65p appears as 



K 2 = 



U2k B T 2 



Mf 2 (47r) 



72-1 



xG" 



92 +72 - 1 



-|2/(3<J2-2) 



(3n - 2) q 



(66) 



As in clusters of galaxies the electron number density is 
typically 10~ 2 ~ 10~ 4 cm -3 (e.g., Fabian 1994; Cavaliere 
& Fusco-Femiano 1978; Nulsen et al. 2005b) and the gas 
mainly consists of protons, electrons, and a particles (nuclei 
of helium atoms), the particle number density is also of this 
magnitude. The mean molecular weight for ICM in galaxy 
clusters is about 0.6 g/mol (e.g., Cavaliere & Fusco-Femiano 
1978) and the mass of hot ICM in galaxy clusters is of the 
order of 1O 13 M (e.g., Peres et al. 1998). Typically, the ICM 
temperature T2 varies in the range of ~ 10 7 — 10 8 K (e.g., 
Fabian 1994). For an observational input of these different 
parameters, we can then estimate the typical range of K2 



5 As an example, we take K2 = ^Q,2 in this subsection for the 
hot ICM on the downstream side of an ICM shock. 



Table 1. Values of K2 in SI unit for different sets of model pa- 
rameters. We take the mean molecular weight of the typical ICM 
in galaxy clusters to be 0.59 g/mol, corresponding to the total 
mass ratio of protons to a particles being 3. Then the ICM mass 
density P2 can be calculated from N e , the electron number den- 
sity in the hot ICM, and K2 is calculated from equation Il66|l . The 
temperature T2 is in unit of keV and the electron number density 
N e is in unit of cm -3 . The enclosed mass M2 for ICM is in unit 
of 10 13 M Q (Mq is the solar mass). All values of K 2 in Table 1 
are in SI unit. 



n 


72 


T 2 




M 2 


K 2 


0.8 


1.31 


6 


0.004 


1 


1.11 x 10 21 


0.8 


1.31 


7 


0.004 


1 


2.68 x 10 21 


0.8 


1.31 


7 


0.008 


1 


7.84 x 10 20 


1 


1.405 


8 


0.002 


1 


3.72 x 10 s 


1 


1.405 


7 


0.002 


1 


6.93 x 10 s 


1 


1.405 


8 


0.001 


2 


1.37 x 10 9 


1.07 


1.42 


5 


0.001 


1 


9.29 x 10 6 


1.07 


1.42 


6 


0.001 


1 


3.98 x 10 6 


1.07 


1.42 


6 


0.003 


1 


3.41 x 10 7 



values (in SI unit). To be specific, we estimate relevant pa- 
rameters for our two-fluid model calculations and the results 
are summarized in Table [1] above. 

3.2 The ICM Temperature Profile 

According to the ideal gas law P2 = piksTi/ 'po, where po 
is the mean molecular mass, fcs is the Boltzmann constant 
and p2 is the ICM mass density, the ICM temperature T2 in 
our self-similar flow model is given by 

T 2 (r,t) = ^£ t 2n - 2 K 2q2 a q 2 2+12 - 1 x 2q2 (Knx~V2) q2 . (67) 
Kb 

For x — > + , the asymptotic behaviour of ICM temperature 
is then 

MO-^2 £2n-2 K 3 92j 492+72-l n <72 a .2-2/n , 
kn 2 



T 2 



(68) 



which diverges as r — + + for n < 1. As the radial flow speed 
Vi(x) is small compared with x as x — > +00 (see equations 
{30} , (J3TJ), (fDTTj) . (DT2)) and the reduced mass density on (x) 
has similar asymptotic scaling behaviours at both large and 
small x (see equations (J2EJ), J2U) and ipgjl). the ICM tem- 
perature T2 has the similar asymptotic scaling behaviours at 
both large and small x because T2 is closely related to ther- 
modynamic variables P2 and p2 as well as the enclosed mass 
M2. Therefore for n < 1, the ICM temperature T2 decreases 
with increasing radius r; in the limit of n — > (2/3) + , the 
limiting temperature scaling would be ~ r _1 . In fact, qual- 
itatively similar temperature profiles have been observed in 
several galaxy clusters (e.g., galaxy clusters A2256, A2319, 
A665 in Markevitch 1996; and galaxy cluster A520 in Marke- 
vitch et al. 2005). Both temperature and mass density di- 
verge as x — > + for n < 1; we need to introduce a sensible 
reference radius to cutoff around the cluster centre. 

With n = 1, expression (I68|l for ICM temperature T2 
remains constant for a certain radial distance around the 
centre. This kind of solution actually represents galaxy clus- 
ters with roughly constant temperature. This kind of galaxy 
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clusters has also in fact been observed, such as galaxy cluster 
A2204 in Peres et al. (1998) and galaxy clusters A2199 and 
2A 0335+096 in Sanders & Fabian (2006). For the n = 1 
case of our two-fluid model, solutions are not 'isothermal' in 
general because 74 ^ 1; and this clearly differs from previ- 
ous isothermal analyses with n = 1 and 7=1 (e.g., Tsai & 
Hsu 1995; Lou & Shen 2004; Lou 2005; Bian & Lou 2005). 
Of course, we could set 71 = 1 and 72 7^ 1 for an 'isother- 
mal' dark matter flow and a nonisothermal gas or vice versa 
or 71 = 72 = 1 for two 'isothermal' flows and so forth. In 
particular, we emphasize that this n = 1 case fundamen- 
tally differs from the condition n — 1 in Suto & Silk (1988), 
because their specific entropy is not conserved along stream- 
lines with n — 1. 

For n > 1, the ICM temperature T2 increases with in- 
creasing radius as indicated by expressions (|67|l and (I68p . 
Temperature profiles in many clusters of galaxies do appear 
to behave in this manner qualitatively in a certain radial 
range, such as galaxy clusters Ms0735.6 + 7421 (e.g., Mc- 
Namara et al. 2005), Perseus (e.g., Sanders & Fabian 2007) 
and A2052 (e.g., Blanton et al. 2001). 



3.3 Energetics of the Coupled Two-Fluid System 

Parallel to the single flow system of an isothermal gas (e.g., 
Tsai & Hsu 1995), the energy of our coupled polytropic two- 
fluid system consists of three parts, namely, the gravitational 
potential energy denoted by E grav , the kinetic energy of 
two fluids denoted by Et,i and the thermal energy of two 
fluids denoted E t h,i, where the subscript i = 1, 2 refer to 
fluid 1 (dark matter) and fluid 2 (hot ICM) respectively. 
In conventional scenarios, AGN activities are sustained by 
accretions of baryon matter onto a supermassive black hole 
(SMBH). Conceptually, it is also physically sensible to think 
of accretions of both baryon matter and dark matter onto 
a SMBH for AGN activities (e.g., Hu et al. 2006), although 
only radiations from ICM can directly reach us. Therefore 
in addition to radiative losses from electrons we observe, the 
outburst energy of an AGN should also involve the energies 
associated with hot gas and dark matter. The energy within 
the shock radius contains the outburst energy of AGN and 
the original energy in the system prior to an AGN. We may 
calculate the energy within the shock radius to estimate the 
order of the outburst energy released by an AGN. 

The gravitational energy of our coupled two-fluid sys- 
tem within a radial range between r m i n and r ma x (the re- 
spective dimensionless self-similar variables are x m i n and 
Xmax for a specified time t) is simply given by 



E 



F 



G(Mi + M 2 



(pi + /52)47rr dr 



(69) 



With self-similar transformation ([5]), we then derive 

5 7^5/2, 5n-4 



Egrav 



p x 

X ., 



(3n - 2)G 





V 2 \ 


^nx - 






K J. 



{a 1 +a 2 )dx , (70) 



where x = x\ for dark matter halo. The kinetic energy for 
fluid i within the radial range between r minti and r max ,i (the 
respective dimensionless self-similar variables are x m i n ,i and 



x m ax,i for a specified time t) is simply given by 

,2 



Eh,i 



piU 



- 4-7rr dr 



2G 



2 2 , 

oaVi Xi dXi 



(71) 



The thermal energy of fluid i within the same radial range 
of r min ,i and r max ,i is simply given by 



EthA = 



Jr mini (7i-l) 



47rr dr 



1 / x 2 a li m qi dx- 



(72) 



As the polytropic sound speed ai in fluid i is 

ai=[ C f i -) =( 7 i-J K't , (73) 

the common dimensional coefficient in integrals (170|l — (1721) 
can be expressed in terms of the sound speed as 



f (-1) 



-5/2 



(74) 



^5/2 i 5n_4 

G 

In our model, the sound speed is not a constant, which 
differs from the isothermal model of Tsai & Hsu (1995) for 
a single gas. 

The total energy of the two-fluid system within the ra- 
dial range between r m i„ and r ma x is then 



Etotal = E gr av + Ek,i + Eth,: 



(75) 



At a certain reference time t\ , we can calculate the total 
energy E to tai,o within a radial range between r min and r max , 
where neither shock arrives. After a certain time lapse t 2 
when both shocks have passed through the radial range un- 
der consideration, we can calculate the total energy E to taij 
within the same radial domain. Then the energy difference 
Etotai, f — Etotai.o is the energy input from shock flow. Fur- 
thermore, the mean input power during this time interval 
t 2 — ti can be estimated by 



Vtotal = (Etotal, f — Etotal,o)/(t 2 — tl) 



(76) 



In addition to radiation losses from the central region, this 
mean input power can be used to estimate the power from 
AGN activities. In our model framework, the total energy in- 
cludes several parts. The X-ray luminosity of galaxy clusters 
inferred by observations (e.g., Nulsen et al. 2005a; Gizani & 
Leahy 2004) should come from the thermal energy part of 
hot ICM. In the estimated outburst energies of AGNs with 
certain models (e.g, Nulsen et al. 2005a; Gizani & Leahy 
2004), only the fraction of energy transmitted to the ther- 
mal reservoir of hot ICM is considered. In our scenario, the 
part of released energy from AGNs is also transmitted to the 
dark matter as well as the gravitational and kinetic energies 
of hot ICM. Therefore, our mean power and energy are more 
than just the thermal energy input to hot ICM. 

3.4 Solution Examples of 2/3 < n < 1 

In this case of 2/3 < n < 1, features of our global semi- 
complete solutions are summarized as follows. When x — > 
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+00, Vi(x) goes to zero with a a.i(x) scaling of x 2// ". In the 
other limit of x — > + , Vi(x) goes to zero with a on(x) scaling 



of x 



-2/r. 



At large x, both fluids may have various combi- 



nations of asymptotic behaviours of outflow, breeze, static 
envelope, contraction, inflow as the two shock positions (or 
outgoing speeds) vary. Through numerical exploration, we 
note that as the shock position decreases from the vicinity 
of the singular surface to smaller x in fluid i, the correspond- 
ing outflow speed Vi(x) at large x tends to decrease. If the 
shock position is reduced further, the fluid at large x may 
become static and then turn into a contraction or an inflow. 
As implied by the limit of our parameter regime for scaling 
index n, the most steep density scaling both at large and 
small x would be ~ x~ 3 (Wang & Lou 2007). This limiting 
density scaling can be systematically tested against obser- 
vations of galaxy clusters. 

Here, we offer a few solution examples with the set of 
parameters {n, 71, 72, ft} = {0.8, 1.3, 1.31, 0.02}. The 
two coefficients A\ and A 2 can be readily calculated from 
equation (|18|) and the value of index k is determined by 
quartic equation (|26[) ( we choose the root k > 1 for the 
consistency of series expansion analysis) ; their specific values 
are {Ai, A2, k} — {1.76, 0.35, 1.99} correspondingly. For 
different initial values and shock positions, the global semi- 
complete solutions for radial flow speeds Vi(x) (scaled for a 
clear presentation) are shown in Figures [T] and [5] numerical 
values of the relevant solution parameters are listed in Table 
[2] In short, for these four numerical solutions 1, 2, 3 and 4, 
the nine parameters to determine a dimensionless solution 
are: {n, 71, 72, ft, Li, L 2 , x ini , x d ,i, x d ,2}={0-8, 1.3, 1.31, 
0.02, 0.06, -0.00638, 6 x 10~ 9 , 342, 0.55} for solution 1 
{n, 71, 72, ft, Li, L 2 , Xmi, Xd,i, x d}2 }={0.8, 1.3, 1.31 
0.02, 0.06, -0.00638, 6 x 10" 9 , 1.8, 0.63} for solution 2 
{n, 71, 72, ft, Li, Z/2, Xini, x d ,i, x dy2 }={0.8, 1.3, 1.31 
0.02, -0.673, 0.0716, 6.8 x 10~ 10 , 0.35, 50} for solution 3 
{n, 71, 72, ft, L\, L 2 , x ini , x d ,i, a;d,2}={0.8, 1.3, 1.31 
0.02, -0.0051, 0.000542, 1.6 x 10~ 8 , 0.8, 5} for solution 4 
It should be noted that parameter Xi„i is not intrinsic to the 
physical description but must be properly chosen within a 
certain range for reliable numerical solutions matched with 
quasi-static asymptotic solutions at small x. 

Solutions such as example solution 2 catch our special 
attention. Such solutions demonstrate a possible scenario 
that when x — > +00, the hot ICM represented by fluid 2 can 
flow outward while the dark matter halo represented by fluid 
1 gradually contracts. As the dark matter halo dominates 
in clusters of galaxies in terms of mass (the enclosed mass 
ratio of dark matter halo to hot ICM varies from 4 to 10 in 
typical clusters of galaxies, e.g., Peres et al. 1998), one can 
readily understand the behaviour of this kind of dynamic 
solutions. Since the hot gas in this solution can overcome 
the overall gravity (dark matter and hot ICM together) and 
flow outward at large distance, we naturally refer to such gas 
outflow as the galaxy cluster wind. The driving energy of the 
galaxy cluster wind may come from the gravitational core 
collapses or AGN activities or merging of galaxies around 
the centre of galaxy clusters; a self-similar phase gradually 
emerges after a while of dynamic evolution. For x — » + , the 
flow velocity of dark matter halo is outward. The outflow of 
dark matter may be propelled by violent outbursts of energy 
around the cluster centre occurred at earlier times. 

In fact, hydrodynamic models for clusters of galaxies 



Table 2. Parameters for numerical solutions 1, 2, 3 and 4 pre- 
sented in Figures[T]and[2]with 2/3 < n < 1. In the first column on 
the left, 'No.' is the numeral label to distinguish different example 
solutions 1, 2, 3 and 4. 'Type 1' is the type of flow solutions of 
fluid 1 when x — > +00, and 'Type 2' is the type of flow solutions 
of fluid 2 when x — > +00. Here, x d 1 and x d 2 are the independent 
self-similar variables on the downstream sides of shock positions 
for fluids 1 and 2, respectively. Parameters E\ and E 2 are re- 
spectively the coefficients of a\(x) and ce 2 (x) at a large enough as- 
according to asymptotic solutions J28D and 11291 . Parameters Hi 
and H 2 are respectively the velocities v± (x) of fluid 1 and v 2 (x) of 
fluid 2 at a large enough x according to asymptotic solutions J30H 
and 13H . Parameter Xi n i is the initial value of x for a numerical 
integration. As this table is too long horizontalwise, we break the 
table in two parts and stack them together. 



No. 


Type 1 


Type 2 


Lx 


L 2 




1 


outflow 


outflow 


0.06 


-0.00638 


6 x 10~ 9 


2 


inflow 


outflow 


0.06 


-0.00638 


6 x 10" 9 


3 


outflow 


inflow 


-0.673 


0.0716 


6.8 x 10" 10 


4 


inflow 


inflow 


-0.0051 


0.000542 


1.6 x 10~ 8 




%d,2 


E 1 


E 2 


Hi 


H 2 


342 


0.55 


6.04 x 10 s 


0.432 


1471 


178 


1.8 


0.63 


2.2 


0.21 


-0.13 


0.034 


0.35 


50 


6.22 


877 


48.5 


-7.61 


0.8 


5 


6.99 


2.67 


-3.44 


-0.0269 



solution 1 and solution 4 
solution 2 



solution 3 



Figure 1. Scaled radial flow velocity of fluid 1 (dark matter halo) 
at small x with 2/3 < n < 1. The relevant model parameters are 
summarized in Table [2] for solutions 1, 2, 3 and 4. Note that for 
a certain range of initial distance, vi(x)/x k remains a constant, 
showing that these solutions for v\(x) are indeed quasi-static for 
small x. For all these solutions 1, 2, 3 and 4, we have the same 
{n, 71, 72, ft} = {0.8, 1.3, 1.31, 0.02}. Solutions 1 and 2 (heavy 
solid line) coincide in the regime of small x. The corresponding 
solution behaviours of v 2 (x)/x k (for the hot ICM) are shown in 
Fig.H 
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Figure 2. Scaled radial flow velocities of fluid 2 (hot ICM) at 
small x with 2/3 < n < 1; the relevant model parameters are 
listed in Table [2] Note that for a certain range of initial distance, 
V2(x)/x k remains constant, showing that these solutions 1, 2, 
3 and 4 for V2(x) are indeed quasi-static for small x. Together 
with Figure [l] for fluid 1 (dark matter halo), we have succeeded 
in constructing solutions with both fluids being quasi-static at 
small x. For all these solutions 1, 2, 3 and 4, we have the same 
{n, 7i, 72, k} = {0.8, 1.3, 1.31, 0.02}. Solutions 1 and 2 (heavy 
solid line) coincide in the regime of small x. 



Figure 3. Negative radial flow velocities of hot ICM (solid curve) 
and dark matter halo (dashed curve) of solution 2 with t = 7.59 X 
10 15 s and Ki = 5.50 X 10 20 SI unit on the downstream side of a 
shock in hot ICM. The shock position in the hot ICM is at ~ 48.38 
kpc while the shock position of dark matter is at ~ 138.23 kpc. 
The initial parameters for the numerical integration are listed in 
Table [2] (see those for solution 2). The corresponding radial flow 
velocities at small x of this solution are shown in Figs, [l] and [2] 
(solution 2). 



have been proposed earlier. In Gunn & Gott (1972), a theory 
of infall of materials into clusters of galaxies was developed 
and they applied it to the growth of galaxy clusters and 
the generation of hot ICM. Bertschinger (1989) discussed 
the time-dependent evolution of cooling flows in clusters of 
galaxies and emphasized that although the size of a cooling 
flow region will increase with time, the gas material itself 
does not go outward; this mechanism is referred to as 'cool- 
ing waves'. In particular, steady winds from a galaxy clus- 
ters (e.g., the Coma cluster) have been discussed by Yahil 
& Ostriker (1973). Their result indicates a mass loss rate of 
~ 10 3 — 1O 4 M0 per year. Their model considered a steady- 
state wind without involving shocks. In our dynamic model, 
we can construct self-similar solutions such as solution 2 to 
explore the dynamics of galaxy cluster winds with shocks 
and possible physical consequences. 

For quasi-static solution 2 in small x shown in Figures 
1 and 2, we adopt estimates for the physical parameters 
* = 7.59 x 10 15 s and K 2 = 5.50 x 10 20 SI unit on the 
downstream side of the ICM shock. In this example, the 
shock position in the hot gas is r s ,2 = 48.38 kpc (shock 
positions vary from tens of kpc to hundreds of kpc in clusters 
of galaxies and a similar shock position has been observed 
in the Perseus cluster; e.g., Fabian et al. 2006). We only 
apply our model to radii less than several Mpcs, say ~ 3 
Mpc (e.g., the size scale of the Hydra A cluster is ~ 3 Mpc; 
see e.g., Taylor et al. 1990). The corresponding flow velocity 
profiles are displayed in Figure [3] and the profile of electron 
number density N e is displayed in Figure [4] In order to get 
the electron number density, we take the mean molecular 
weight to be ~ 0.59 g/mol (see Cavaliere & Fusco-Femiano 
1978). In this example, the radial outflow velocity of hot gas 
at 3Mpc is ~ 202 km s _1 while the radial inflow velocity of 
dark matter halo at 3Mpc is ~ 20.3 km s _1 . The analytical 
expression for the polytropic sound speed ratio can be found 




10' io z 10 3 



r (kpc) 

Figure 4. Electron number density profile for t = 7.59 x 10 15 s 
and K 2 = 5.50 X 10 20 SI unit on the downstream side of a shock 
for solution 2 (see Figs. 1 and 2 for details of the quasi-static 
solution behaviour at small x). The corresponding global radial 
flow velocities of this solution are shown in Fig. [3] Here we take 
the typical mean molecular weight in clusters of galaxies to be 
0.59 g/mol. As the electron number density can be obtained by 
high-resolution X-ray imaging observations, we give the model 
electron number density so that comparisons and tests can be 
made. 

in Appendix G and the sound speed ratio in this specific 
example is displayed in Fig. [6] 

According to expression (|67[) . the temperature of hot 
ICM on the downstream side of the shock is ~ 0.85keV 
while the temperature on the upstream side is ~ 0.49 keV. 
The enclosed mass of hot gas at 0.5Mpc is about M2 = 
2.4 x 10 12 Mq; the enclosed mass ratio of dark matter halo 
to hot gas at 0.25Mpc is ~ 7.58 and at 0.5Mpc is ~ 7.47. 
By expression Q440 . the outgoing shock speed in the hot 
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Figure 5. The enclosed mass ratio M2/M1 as a function of r at 
different times of solution 2. The abscissa is the radius in the unit 
of virial radius r2oo (defined in the paragraph before subsection 
3.1) while the ordinate is the enclosed mass ratio M2/M1 between 
the hot ICM and dark matter halo. The three curves are the 
enclosed mass ratio of solution 2 with the same parameter K2 = 
5.50 X 10 20 in SI unit on the downstream side of a shock but with 
different times as indicated along the curves. As for different times 
the virial radius T200 is a little different, the same position in the 
space corresponds to different radius in the figure in principle. 
However, the differences between the virial radius at different 
times are so small that we can treat one specific point of the x axis 
as the same position in the space at different times approximately. 
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Figure 6. The radial profile for the polytropic sound speed ratio 
a 2 /a± of fluid 1 (dark matter) to fluid 2 (hot ICM) at time t = 
7.59 X 10 15 for quasi-static solution 2 shown in Figs. 1 — 5. The 
first jump on the left at about r = 48.38 kpc is due to the shock 
in fluid 2 (hot ICM) and the second jump on the right at about 
r = 138.23 is due to the shock in fluid 1 (dark matter halo). When 
x — » + , this ratio is about a\/a,2 = 0.99 and when x — » +oo this 
ratio is about a\ja2 = 1.216. 



ICM is u St 2 = nr at 2/t = 157.6 km s . For a shock located 
at a radius of 48.38 kpc with a reference timescale of t — 
2.4 x 10 8 yr, the shock position by equation (|45[) is r s ,2 = 
9.54 x 10 -6 i ' 8 where r Sj 2 is in the unit of kpc and time 
t is in the unit of year. For example, during a timescale of 
~ 10 10 yr, this shock would travel to a radius of ~ IMpc 
within a cluster of galaxies. The enclosed mass ratio M%/Mi 
between the hot ICM and the dark matter of this solution 
at different radii is shown in Fig. (the solid line with time 
t — 2.4 x 10 8 yr). If we adjust the timescale, the solution 
will evolve in the self-similar manner. We choose two other 
timescales as examples with other parameters of solution 2 
fixed and the result is also shown in Fig. [5] Note that with 
increasing time, the minimum of this ratio moves towards 
larger radii which is the result of accretion of hot gas and 
outflow of dark matter around the centre. The ratio remains 
nearly independent of time at the centre and at large radii. 

As an example of conceptual exercise, we calculate at 
the reference time fi = 7.59 x 10 15 s the energies within the 
radial range ~ 138 — 150.8 kpc undisturbed by shocks. The 
gravitational potential energy is E gra v,o = —2.66 x 10 60 erg 
by expression (|70p . the kinetic energy of fluid 1 is Ek,i,o = 
4.94 x 10 57 erg by expression (|7ip. the kinetic energy of 
fluid 2 is E k ,2, = 4.26 x 10 56 erg by expression (fTl), the 
thermal energy of fluid 1 is E t h,i,o = 3.43 x 10 60 erg by 
expression (|72p , and the thermal energy of fluid 2 is E t h,2,o = 
2.72 x 10 59 erg by expression (|72p . According to equation 
([751) . the total energy is then E to tai,o = 1-05 x 10 60 erg. As 
r s ,2 = 9.54 x 10~ 6 t 0,8 , the two shocks have passed through 
the radius 150.8 kpc by the time t 2 = 3.14 x 10 16 s. Then 
at this time t2, the energies within the same radial range 
are: the gravitational energy E grav j — —2.55 x 10 60 erg 
by expression ((70}, the kinetic energy of fluid 1 Et,i,f = 
1.0 x 10 56 erg by expression (|7ip , the kinetic energy of fluid 
2 E ka ,f = 8.32 x 10 56 erg by expression J7TJ, the thermal 
energy of fluid 1 E t h,i,f = 3.41 x 10 60 erg by expression 
(f72)l . and the thermal energy of fluid 2 E t h,2j = 3.09 x 10 59 
erg by expression (I72p . According to equation (|T5f> . the total 
energy is then Etotaij = 1T6 x 10 60 erg. Then according to 
equation ([76} , the mean power of shock flow in this example 
is Vtotai = 4.94 x TO 42 erg s _1 . 

Around a distance of r = 136 kpc and with increasing 
r, the radial flow velocity of hot gas changes from inflow 
to outflow, where we can calculate the total mass accretion 
rate. We take a position r a = 135.98 kpc and the mass ac- 
cretion rate there is then M a ,2 = 4-7rp2W2i~a ~ 27. 1M© yr -1 , 
which is comparable to the mass accretion rates inferred for 
galaxy clusters A3158(P) and A262(P). There is a summary 
of these parameters for different clusters of galaxies, which 
are grossly consistent with our numerical example illustrated 
here (e.g., Peres et al 1998.) 

At a radial distance of ~ 3Mpc, the outflow mass 
per year of hot ICM is ~ M 2 ,o = 4np2U2r 2 = 211M 
yr -1 . Thus the approximate total outflow mass of hot gas 
within a timescale of ~ 2.4 x 10 yr is 5.1 x 10 10 M Q . 
Meanwhile, the inflow mass per year of dark matter is 
~ Mi,o = 4-rvpimr 2 = 140M Q yr -1 . Then the total in- 
flow mass of dark matter within a timescale of 2.4 x 10 8 
yr is approximately ~ 3.36 x 1O 1O M . Here, we take the 
timescale of galaxy cluster winds to be ~ 6 x 10 9 yr. Then 
during this time, we assume for simplicity that there is an 
AGN in a cluster of galaxies every 2.4 x 10 8 yr and the 
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hot gas mass of ~ 5.1 x 10 Mq is lost and the dark mat- 
ter mass of ~ 3.36 x 1O 1O M0 is deposited each time. (That 
means there have been roughly 25 times of AGN activities 
in galaxy clusters since initial galaxy cluster winds formed.) 
Furthermore, when galaxy cluster winds began to form, we 
may presume that the baryon mass fraction is roughly the 
average cosmic baryon fraction in the expanding universe. 
In this example, the total mass of hot ICM within 3 Mpc is 
6.1 x 10 12 A/q and the total mass of dark matter within this 
same radial scale is 4.1 x 1O 13 M0 at present. Then the value 
of baryon fraction fb at present is 0.129 while the value of 
baryon fraction /j, when this cluster of galaxies begins to 
blow a galaxy cluster wind is 0.154. Therefore the baryon 
fraction /(, decreases 16.4% due to galaxy cluster winds and 
the inflow or contraction of dark matter halo in this example 
of illustration. In this scenario, the 'missing baryons' should 
be found in the periphery of galaxy clusters on much larger 
scales. They cool down to lower temperatures and may not 
be easily detected. 

3.5 Solution Examples of n — 1 

When n — 1, the radial flow velocities of both fluids ap- 
proach constant values as x — > +oo; for global semi-complete 
solutions, these values can be negative, zero or positive at 
large x, depending on the choice of relevant parameters and 
positions of shocks. With n = 1, the asymptotic constant 
speeds are Hi and H2 in asymptotic velocity solutions (|30p 
and (|31[) respectively. Physically, these correspond to in- 
flows, breezes or contractions, and outflows, respectively. As 
x — *• + , flow solutions can be matched with the asymp- 
totic quasi-static solution for both fluids. For n = 1, the 
mass density profile is pi — » r~ 2 either as r — » + or as 
r — > +00. We now provide a few examples with parame- 
ters {n, 71, 72, k} = {1, 1.3, 1.405, 4} and correspond- 
ing values of {Ai, A 2 , k} being {267.2, 38.4, 1.89}. Sim- 
ilar to the case of2/3<n<las described in the pre- 
vious subsection, we start the numerical integration with 
a small enough Xini and go to a large enough x max . Once 
vi(x) and V2(x) become constant after numerical integra- 
tions in a certain small x range, we would regard the x max 
as sufficiently large. The results are shown in Fig. [7] to 
Fig. \W\ and the relevant model parameters are summarized 
in Table [3] Specifically for these solutions 5, 6 and 7, the 
nine parameters to determine a dimensionless solution are: 
{n, 71, 72, K, Li, L 2 , x lni , x d: i, x da } = {1, 1.3, 1.405, 4, 

2.66 x 10 -3 , -0.0913, 1 x 10~ 8 , 25.1, 6.5} for solution 5; 
{n, 71, 72, k, Li, L 2 , Xini, x d ,i, Xd,2}={l, 1-3, 1.405, 4, 

9.67 xl0~ 4 , -0.0331, 3 xl0" 8 , 15, 9.5} for solution 6; and 
{n, 71, 72, k, Li, L 2 , Xini, x d ,i, Xd,2,}={l, 1-3, 1.405, 4, 
9.67 x 10" 4 , -0.0331, 3 x 10 -8 , 12, 9.5} for solution 7. 

Again, we use solution 6 to further illustrate the process 
of a galaxy cluster wind. Here, fluid 1 represents the dark 
matter halo and fluid 2 represents the hot ICM. We choose 
the parameters t and K2 (the value on the downstream side 
of an ICM shock) as t = 6.32 x 10 15 s and K 2 = 7.34 x 10 s 
SI unit. In this example, the shock positions in the hot ICM 
and in dark matter halo are r s ,2 = 210.75 kpc and 332.76 
kpc, respectively. In galaxy cluster Hydra A, the observed 
shock is located at ~ 211 kpc (e.g., Nulsen et al. 2005b) 
and the shock in galaxy cluster Ms0735. 6+7421 is at ~ 240 
kpc (e.g., McNamara et al. 2005). For this example, we ap- 



Table 3. Model parameters for solutions in Fig. [7]to Fig. HOI with 
n = 1. Here, No. on the left gives the numeric label to distinguish 
different solutions. Type 1 is the type of solutions for fluid 1 (dark 
matter halo) when x — > +00 and Type 2 is the type of solutions 
for fluid 2 (hot ICM) when x — > +00. Here, Xd 1 and Xd 2 are 
the independent self-similar variables on the downstream sides 
of the respective shock positions for fluids 1 and 2. Ei is the 
coefficient of ai and E2 is the coefficient of a.2 at a sufficiently 
large x (see asymptotic solution 1281 and I29H . Hi is the velocity of 
fluid 1 at a sufficiently large x and H2 is the velocity of fluid 2 
at a sufficiently large x (see equations 1301 and I3U . As this table 
is too long horizontally, we break this table into two parts and 
stack them together. 



No. Type 1 Type 2 L x L 2 x ini 

5 outflow outflow 2.66 x 10" 3 -0.0913 1 x lO - ' 

6 inflow outflow 9.67 x 10~ 4 -0.0331 3 x 10~ 

7 inflow inflow 9.67 x 10~ 4 -0.0331 3 x 10~ 



Xd,l x d,2 Ei E 2 Hi H 2 

25.1 6.5 487 33.3 7.84 49.5 

15 9.5 265 30.5 -0.811 10.34 

12 9.5 200 28.2 -5.00 -2.34 




solution 5 



15 solution 6 and 

solution 7 

, _/_._. 



'l<T* 10"" 10"" 10~ ! 10° 10 ! 

x 

Figure 7. Scaled radial flow velocities of fluid 1 (dark matter) at 
small x when n = 1. The relevant model parameters are summa- 
rized in Table [3] Note that for a certain range of radial distance 
around small x, vi/x k remains constant, showing that these solu- 
tions 5, 6 and 7 for vi are indeed quasi-static solutions as derived 
analytically. For all these three solutions 5, 6 and 7, we have the 
same parameter set {n, 71, 72, k} = {1, 1.3, 1.405, 4}. The 
corresponding flow velocities of fluid 2 (hot ICM) are shown in 
Fig. E] 



ply our model to the radial domain of ~ 1 Mpc (some rich 
clusters of galaxies may have radial size scales of ~ 1 — 2 
Mpc; e.g., Bahcall 1996). The radial flow velocity profiles 
are shown in Fig. [11] while the profiles of electron number 
density and ICM temperature are displayed in Figs. [12] and 
1131 respectively. Actually the temperature profile similar to 
Fig [13] in the entire radial range has been observed in galaxy 
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Figure 8. Radial flow velocities of fluid 2 (hot gas) at small x 
when n = 1. The relevant model parameters are listed in Table 
\3\ Note that for a certain range of initial distance, V2/x k remains 
constant, which shows that these solutions of V2 are indeed quasi- 
static solutions. Together with Fig. [7] for radial flow velocities of 
fluid 1 (dark matter halo), we have succeeded in constructing so- 
lutions with both fluids being quasi-static at small x. For all these 
solutions, we have the same {n, 71, 72, k} = {1, 1.3, 1.405, 4}. 



Figure 10. Radial flow velocities of fluid 2 at large x when 
n = 1 . The relevant model parameters are listed in Table [3] with 
{n, 71, 72, re} = {1, 1.3, 1.405, 4}. As shown in this figure, i>2 
remains constant for these three solutions at large x, indicating 
that the x max we adopt is sufficiently large for the asymptotic 
behaviour of V2(x). For solutions 6 and 7, the shock positions are 
the same in gas but are different in dark matter halo (see Fig. 
[9j . Together with Fig. [9] we have succeeded in constructing the 
semi-complete solutions with n = 1. The radial flow velocities of 
fluid 2 at small x for these solutions are shown in Fig. [8] 
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Figure 9. Radial flow velocities of fluid 1 (dark matter halo) at 
large x for n = 1. The relevant model parameters are listed in 
Table |3] with {n, 71, 72, k} = {1, 1.3, 1.405, 4}. As shown 
in this figure, vi remains constant for these three solutions at 
large x, indicating that the x max we adopt is large enough for 
the asymptotic behaviour of v\(x). The radial flow velocities of 
fluid 1 at small x for these three solutions are shown in Fig. 
The shock of solution 5 locates at the turning point x = 25.1, 
jumping from v\ = 8.645 to v\ = 8.589 and the shock of solution 
6 locates at x = 15, jumping from v\ = 0.99 to v\ = 0.3. 
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Figure 11. Radial flow velocities of hot ICM (solid curve) and 
dark matter halo (dashed curve) of solution 6 with n = 1, t = 
6.32 X 10 15 and K 2 = 7.34 X 10 8 SI unit in the downstream 
side of a shock in hot ICM. The shock positions in hot ICM and 
dark matter halo are at ~ 210.75 kpc and at ~ 332.76 kpc. The 
initial parameters for numerical integration are listed in Table [3] 
(solution 6). Details of the dimensionless velocities of this solution 
are shown in Figs. [7] \E\ [9] [TO] (solution 6). 



cluster NGC1275 (e.g., Churazov et al. 2003). The ratio of 
polytropic sound speeds (see Appendix G) in two fluids in 
this example is shown in Fig. 1151 In this example, the radial 
outflow velocity of hot ICM at IMpc is ~ 97.4 km s _1 and 
the radial inflow velocity of dark matter halo at 3 Mpc is 
~ 66.76 km s - . The temperature of hot ICM downstream 
of the shock is ~ 9.56 keV and that of the upstream of 
the shock is ~ 7.9 keV. The enclosed mass of hot ICM at 
0.5Mpc is about M2 — 4 x 1O 13 M0; and the enclosed mass 
ratio of dark matter halo to hot ICM at 0.25Mpc is ~ 8.1 



and at 0.5Mpc is ~ 9.3. The travel speed of the outgoing 
shock in the hot ICM is u s 2 — nr a ,2 A = 1.03 x 10 3 km 
s _1 . For a shock located at 210.75 kpc with a timescale of 
t — 2 x 10 8 yr, the shock position (see expression (|45p ) can 
be calculated by r 3j 2 = 1.05 x 10 -6 t where r s .2 is in the 
unit of kpc and time t is in the unit of year. For n = 1, 
the shock speed remains constant. The enclosed mass ratio 
Ma/Mi between the hot ICM and the dark matter of this 
solution 6 shown Figs. 7 — 13 at different radii is displayed 
in Fig. [14] (the solid line with time t = 2 x 10 8 yr). If we 
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Figure 12. The radial profile of electron number density with a 
timescale of ( = 6.32 X 10 15 s and K% = 7.34 X 10 s SI unit on 
the downstream side of a shock in the hot ICM for solution 6 of 
n = 1. Here, we also take the typical mean molecular weight in 
clusters of galaxies to be ~ 0.59 g/mol. The radial flow velocities 
of this solution 6 are shown in Fig. 1111 
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Figure 13. The radial temperature profile of hot ICM with a 
timescale of i = 6.32 X 10 15 s and K 2 = 7.34 X 10 s SI unit on the 
downstream side of a shock in ICM for solution 6 of n = 1. During 
a range of short radial distance from the centre, the temperature 
is a constant of ~ 11 keV and at large x the temperature is ~ 8.2 
keV. The electron number density of this solution 6 is shown in 
Fig. 1121 and the corresponding radial flow velocities are shown in 
Fig. OH 



adjust the timescale t, solution 6 will evolve in a self-similar 
manner. We choose two other timescales t = 10 8 yr (dashed 
curve) and 4 x 10 7 yr (dash dotted curve) as two additional 
examples with other parameters of solution 6 unchanged and 
the results are shown in Fig. 1141 With increasing time t, the 
minimum of this enclosed mass ratio moves towards larger 
radii which is the result of accretion of hot ICM and outflow 
of dark matter around the centre. This enclosed mass ratio 
is nearly independent of time t around the centre and at 
large radii. 

As another example of conceptual exercise, we calculate 
at a reference time t = 6.32 x 10 15 s the energies within radial 
range ~ 332.7 — 502 kpc undisturbed by shocks below. The 




radius/r„. 



Figure 14. The radial profile of enclosed mass ratio M2/M1 at 
different times of solution 6 shown in Figs. 7 — 13. The abscissa 
is the radius in unit of virial radius r2oo and the ordinate is the 
enclosed mass ratio M2/M1 between the hot ICM and dark mat- 
ter. The three curves are enclosed mass ratio of solution 6 with 
the same parameter K2 = 7.34 X 10 8 SI unit on the downstream 
side of an ICM shock and three different times t marked along the 
curves. As for different times the virial radius r200 is a somewhat 
different, the same position in space corresponds to different radii 
in this figure. However, these differences among the virial radii at 
different times are so small that we may treat one specific point 
of the abscissa as the same position in space at different times 
approximately. 




Figure 15. The radial profile for the ratio ai/02 
sound speeds in fluid 1 (dark matter) to fluid 2 (hot 
t = 6.32 X 10 15 for quasi-static solution 6 shown in 
The first jump on the left at about r = 210.75 kpc 
shock in fluid 2 and the second jump on the right at 
kpc is due to the shock in fluid 1. When x — » 0^ 
speed ratio is about a\/a,2 = 0.962 and when x — » 4 
is about ai/a 2 = 1.11. 



of polytropic 
ICM) at time 
Figs. 7 - 14. 
is due to the 
about 332.76 
", this sound 
00, this ratio 
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gravitational energy is E grav , — — 9.74x 10 erg (see equa- 
tionlTO)!. the kinetic energy of fluid 1 is E k ,i,o = 7.39 x 10 59 
erg (see equation [71}, the kinetic energy of fluid 2 is Ek,2,o = 
8.42 x 10 60 erg (see equation [7T}, the thermal energy of fluid 
1 is E t h,i,o = 1-97 x 10 64 erg (see equation [72} , the thermal 
energy of fluid 2 is E t h,2,o = 1.05 X 10 63 erg (see equation 
[72}. According to equation (|75[) . the total energy is then 
given by E tota l,o = 1.11 X 10 64 erg. As r s , 2 = 1.05 x 10~ 6 t, 
the two shocks have passed through the radius 502 kpc by 
time t2 — 1.50 x 10 16 s. Then at this time the energies 
within the same radial range are: the gravitational energy 
Egravj = —1-01 x 10 64 erg (see equation[7D}, the kinetic en- 
ergy of fluid 1 Ek,ij = 5.78 x 10 58 erg (see equation [71} . the 
kinetic energy of fluid 2 Ek,2,f = 7.76 x 10 59 erg (see equa- 
tion|7T}, the thermal energy of fluid 1 E thx} = 1.97 x 10 64 
erg (see equation[72}, the thermal energy of fluid 2 Eth,2,f = 
1.73 x 10 63 erg (equation [72}. According to equation ([75} . 
the total energy is then E to taij = 1-13 x 10 64 erg. Then 
according to equation (|76[) . the mean power of shock flow in 
this illustration example is Vtotal = 2.94 x 10 46 erg s _1 . 

Around a radial distance of 695.5 kpc and with increas- 
ing r, the flow velocity of hot gas changes from inflow to out- 
flow, where we can calculate the total mass accretion rate. 
We take a radius of r a — 695 kpc and the mass accretion 
rate is then M a ,2 = ^p2U2r 2 a ~ 193M yr -1 , comparable 
to the mass accretion rates in galaxy clusters A85(P) and 
A644(P) (e.g., Peres et al. 1998). 

At a radius of ~ IMpc, the outflow mass per year of hot 
ICM is ~ M 2 ,o = 4np 2 u 2 r 2 = 8.11 x lO 3 A/ yr~\ which is 
comparable to the steady wind result of Yahil & Ostriker 
(1973). Thus the approximate total outflow mass of hot gas 
within a timescale of ~ 2 x 10 8 year is ~ 1.6 x 1O 12 M . 
Meanwhile, the inflow mass per year of dark matter is 
~ Mi, = iivpiUir 2 = 4.94 x 10 4 A/ Q yr -1 . Then the to- 
tal inflow mass of dark matter within a timescale of 2 x 10 8 
yr is approximately 1O 13 M . Here, we take the timescale of 
galaxy cluster winds to be 3 x 10 9 yr and the AGN activ- 
ity occurs every 2 x 10 8 yr. (That means there have been 
roughly 15 times of AGN activities in galaxy clusters since 
the initial wind formed.) In this example, the total mass of 
hot ICM within 1 Mpc is 8.1 x lO 13 A/ and the total mass of 
dark matter within this same radial scale is ~ 7.4 x 1O 14 M 
at present. Now consider the mass loss of hot ICM due to 
galaxy cluster winds and inflow of dark matter. The total 
mass of hot ICM within this radial domain when the galaxy 
cluster winds began to form was 1.05 x 1O 14 M and the total 
mass of dark matter within this same radial domain when 
the galaxy cluster winds began to form was 5.96 x 10 14 M Q . 
Therefore the baryon fraction /& at present is 0.098 and its 
value when the cluster began to form winds is 0.15. The 
value of fb decreases 34.7% due to galaxy cluster winds and 
inflow of dark matter halo in this example of n — 1. 

3.6 Solution Examples of n > 1 

With n > 1 in our model, radial flow velocities of hot gas 
and dark matter will tend to diverge for Hi and H2 7^ 
as we integrate towards x — » +00 (see asymptotic solutions 
([30} and ([31} towards large a;) with quasi-static asymptotic 
solutions around the cluster centre. This case of n > 1 can 
be relevant as the density profile scales as p oc r~ 2 ^ n while 
in many clusters of galaxies the electron number densities 
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Figure 16. The scaled radial flow velocities as x — > + for conver- 
gent solutions at large x when n > 1 and with Hi = H2 = 0. The 
model parameters at small x are {n, 71, 72, k, Xini, L\, L2] = 
{1.07, 1.315, 1.42, 0.1, 1 x 10" 8 , 3.65 x 10" 5 , -2.94 x 10" 5 }. 
This figure shows that the solution approaches quasi-static con- 
figuration for both fluids with this set of parameters at small x. 
These parameters are kept the same as we try to find the conver- 
gent solution. So the solution we find here is still quasi-static at 
small x. 

observed scale as a power law with power indices between 
— 1 and —2 (e.g., clusters of galaxies A2204, A2052 and 
Ms0735. 6+7421). So for a physical system, we only apply 
our model to a finite size of the order of ~ 1 Mpc, which is 
the size of a typical cluster of galaxies. Therefore, we have 
considerable interest in those solutions with two coefficients 
Hi and H2 of the diverging terms being zero. We are now in 
a position to show such an example of solutions with n > 1. 

In our model, we choose relevant parameters to be 
{n, 71, 72, k} = {1.07, 1.315, 1.42, 0.1} with the cor- 
responding values of {Ai, A 2 , k} = {1.53 x 10 5 , 2.27 x 
10 4 , 1.81}. The parameters at small x are {2;™, Li, L2} = 
{1 x 10" 8 , 3.65 x 10~ 5 , -2.94 x 10~ 5 }. The radial flow ve- 
locities at small x are shown in Fig. 1161 confirming that the 
solution we find is indeed quasi-static as x — > + . By numer- 
ical exploration, we can choose proper values of independent 
self-similar variables Xd,\ and Xd,2 on the downstream sides 
of respective shock positions and construct a solution with 
Hi = and H2 = (i.e., finite radial flow velocities at 
large a;). With the relevant parameters chosen for small x, 
we vary the shock positions in fluid 1 and fluid 2, respec- 
tively, and compute the coefficients Hi and I/ 2 . When the 
product Vi(x) x~ 1+1 ^ n approaches constant values for either 
values of i after a certain x max , we would regard x as being 
large enough for asymptotic solutions and then evaluate the 
values of Hi and iif 2 . Empirically, we find that for a fixed 
%d,2 value, if we change the value of Xd,i continuously, Hi 
and H2 lie on a perfect straight line. When we vary Xd,2, 
this Hi versus H2 line will move. For a certain Xd,2, this line 
can move across the zero point in the functional relation 
of Hi versus i/ 2 . When this happens, there exists a pair of 
{xd,i, %d,2\ such that Hi — H2 = 0. This corresponds to 
convergent radial flow velocities for both fluids at large x. 
The result is then shown in Fig. [17] and the relevant param- 
eters are summarized in Table [4] 
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Table 4. Model parameters for constructing a solution of n > 1 
with Hi = and H 2 = 0. Here, Xd,i and 2^ 2 are the indepen- 
dent self-similar variables on the downstream sides of the shock 
positions for fluid 1 (dark matter) and fluid 2 (hot ICM), re- 
spectively. Hi is the coefficient of vi(x) and H2 is the coeffi- 
cient of v 2 (x) when x becomes sufficiently large. Other model 
parameters at small x are {n, Li, L 2 } = 

{1.07, 1.315, 1.42, 0.1, 1 x 10" 8 , 3.65 x 10" 5 , -2.94 x 10~ 5 }. 
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Figure 18. The radial flow velocities of hot ICM and dark matter 
halo of convergent solution with n > 1, t = 1.1 x 10 16 s, and a 
downstream K 2 = 5.03 x 10 6 SI unit in the hot ICM. The shock 
position in hot ICM is at ~ 298.6 kpc and the shock position in 
the dark matter halo is at ~ 596.2 kpc. The relevant parameters 
for the numerical integration are {n, 71, 72, n, Xi n i, Li, L2] = 
{1.07, 1.315, 1.42, 0.1, 1 x 10" 8 , 3.65 x 10" 5 , -2.94 x 10" 5 }. 
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Figure 19. The electron number density profile of a convergent 
solution of n > 1 with ( = l.lx 10 16 s and K 2 = 5.03 X 10 6 in SI 
unit in the downstream side of the shock in the hot ICM. Here, 
we also take the typical mean molecular weight in galaxy clus- 
ters to be 0.59 g/mol (e.g., Cavaliere & Fusco-Femiano 1978). 
The other parameters are {n, 71, 72, k, Xi n i, Li, L 2 } = 
{1.07, 1.315, 1.42, 0.1, 1 x 10" 8 , 3.65 x 10" 5 , -2.94 x 10" 5 }. 



Figure 17. The numerical search of a global semi-complete so- 
lution of n > 1 with Hi = and H 2 = 0. As an example, the 
dash-dotted line marked by 285 is the line with x^ 2 = 285. As 
Xd 1 increases, the point goes from the bottom left to upper right 
(values of both Hi and H 2 increase) along a straight line. Other 
straight lines are plotted in the same manner with different val- 
ues of x^ 2 as explicitly marked. The solid, heavy dashed, and 
dotted lines are with x^ 2 = 280, 277, 275, respectively. The 
trend of variation is clear by this numerical exploration. The pa- 
rameters of these four straight lines are summarized in Table [4] 
Other parameters at small x are {n, 71, 72, K, Xi n i, L\, L 2 } = 
{1.07, 1.315, 1.42, 0.1, 1 x 10~ 8 , 3.65 x 10~ 5 , -2.94 x 10~ 5 }. 



From Fig. 1171 we know the existence of a solution with 
convergent flow velocities at large x with Xd, 2 = 280 and 
Xd.i = 559. We shall describe the physical process that this 
solution represents presently. The time t is chosen to be 
t — 1.1 x 10 16 s and the downstream K 2 of hot ICM is chosen 
to be K 2 = 5.03 x 10 6 SI unit. Then the shock position in 
hot ICM is at ~ 298.6 kpc and the shock in the dark matter 
halo is at ~ 596.2 kpc. 

We apply this numerical example to within a radial scale 
of ~ 3 Mpc. The radial flow velocity profiles are shown in 
Fig. 1181 If the two coefficients Hi and H 2 are exactly zero, 
the hot gas velocity is inward and vanishes, while the flow 
velocity of dark matter is outward and vanishes as r — > 
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Figure 20. The temperature profile of hot ICM for a convergent 
solution at large x of n > 1 with t = 1.1 X 10 16 s and K2 = 
5.03 X 10 6 in SI unit in the downstream side of the shock in the 
hot ICM. The temperature increases with increasing radius at 
both small and large x. 

+00. The profile of electron number density is shown in Fig. 
1191 and the temperature profile is shown in Fig. 1201 In this 
example, the inflow velocity of hot gas at 3Mpc is about 
231.4 km s _1 while the outflow velocity of dark matter at 
3Mpc is ~ 32.8 km s~ . Across the shock, the temperature 
of hot ICM on the downstream side is ~ 9.5 keV and that 
on the upstream side is ~ 6.17 keV. The enclosed mass of 
hot ICM at 0.5Mpc is about M 2 = 3.33 x 1O 13 M ; and 
the enclosed mass ratio of dark matter to the hot ICM at 
0.25Mpc is ~ 7 and at 0.5Mpc is ~ 9.5. The outgoing shock 
speed in the hot ICM is u Sj 2 = nr s ^/t — 892 km s _ . For 
a shock located at r = 298.6 kpc with a timescale of t = 
3.8 x 10 s yr, the shock position given by expression (|45[1 can 
be followed by r s , 2 = 2.2 x 1CT 7 t im where r s , 2 is in the 
unit of kpc and time t is in the unit of year. 



4 SUMMARY 

In this paper, we formulate a two-fluid model framework 
of spherical symmetry to explore dynamic behaviours of 
hot ICM and dark matter during the evolution of galaxy 
clusters. In this scenario, the hot ICM and dark mat- 
ter halo are approximated as two polytropic 'fluids' and 
are coupled by gravity. For both 'fluids' in general situ- 
ations, specific entropies are conserved along streamlines 
separately and are related to the enclosed masses. Quasi- 
static solutions for both 'polytropic fluids' can be ob- 
tained and are adopted for sufficiently small radii around 
the central core region at a given time. In order to con- 
struct dimensionless quasi-static solutions in the regime of 
small x, we need to specify nine dimensionless parameters 
{n, 71, 72, k, Li, L 2 , Xi„i, x d ,i, 2^,2} where x ini needs to 
be properly chosen. Two additional dimensional parameters 
K2 and t need to be specified for dimensional solutions for a 
physical description. The ICM temperature is taken to be on 
the order of ~ 10 7 — 10 8 K in typical clusters of galaxies and 
the electron number density is taken to be ~ 10 -2 — 1CP 4 
cm -3 in the typical radial range of kpes to Mpcs. It is possi- 
ble to construct different types of flow solutions with shocks 



in both fluids and with various asymptotic scaling features in 
flow speeds, mass densities, enclosed mass and temperatures. 
In particular, we can construct dynamic solutions for galaxy 
cluster winds and discuss the important problem of 'miss- 
ing baryons' during the evolution of galaxy clusters. There 
are several physical hypotheses in our two-fluid model to 
simplify the mathematical analysis. First, we take the dark 
matter halo as a kind of 'fluid' for simplicity and explore this 
alternative theoretical possibility. Secondly, we assume the 
two-fluid system of galaxy clusters to be grossly spherically 
symmetric with a common centre and hope to catch major 
dynamic flow features on large scales. Thirdly, both 'fluids' 
are assumed to be 'polytropic' in the most general sense of 
entropy conservation along streamlines. In other words, the 
specific entropy distribution is not necessarily constant in 
space and time but are allowed to vary in r and t in gen- 
eral. Finally, dark matter interacts with the hot ICM only 
through gravity. 

Global semi-complete solutions can be constructed to 
pass through the two singular surfaces via shocks in both 
fluids. As large-scale ICM shocks have been identified obser- 
vationally, there must be large-scale flows of hot ICM. As the 
dark matter halo and the hot ICM are coupled by gravity, 
there may be dark matter flows and the possibility of shocks 
in dark matter halo. Such dark matter shocks are character- 
ized by drastic density jumps and sharp rises of velocity 
dispersions and may be detected by utilizing gravitational 
lensing effects. In our model framework, outflows of hot ICM 
in galaxy clusters actually form galaxy cluster winds, which 
is a systematic mechanism of reducing the baryon fraction 
ft,. In our model, the self-similar shocks travel outwards in 
both hot ICM and dark matter halo, respectively. The travel 
speeds of these shocks actually change with time for n^l 
and their radial positions at present can vary from tens of 
kpes to a few Mpc. Due to galaxy cluster winds, the baryon 
fraction ft in galaxy clusters can be ~ 15% — 40% lower than 
the average cosmic baryon fraction in the Universe, which 
may account for the problem of 'missing baryons' in clus- 
ters of galaxies. Physically, these 'missing baryons' should 
reside in the periphery of galaxy clusters in the form of 
warm gas as results of unavoidable radiative cooling. Since 
the lower baryon fraction ft is a generic phenomenon in 
clusters of galaxies, we therefore suggest that galaxy cluster 
winds would be common and frequent during the evolution 
of galaxy clusters. 

The main features of our self-similar polytropic solution 
for two fluids coupled by gravity are summarized below. The 
radial profile of mass density at a given time is pi oc r~ 2//n in 
both fluids for either r — > + or r — *• +00 with 2/3 < n < 2 
in general. That shows that the asymptotic mass density 
profiles of hot ICM and dark matter take the same form 
of power laws with the same index — 2/n. Meanwhile, the 
radial flow velocities of both fluids approach zero as x — > + 
in the form of quasi-static asymptotic solution (I24p and (1251) . 
At large x, the asymptotic flow solution is (1281) . (1291) . (|30[) . 
(pil for 2/3 < n < 2. When 2/3 < n < 1 and x -> +co, 
the asymptotic radial flow velocities of both fluids approach 
zero. When n = 1 and x — > +00, the radial flow velocities of 
both fluids approach constant values, which may be positive, 
zero or negative for various combinations of two fluids. When 
n > 1 and x — > +00, the radial flow velocities of both fluids 
become divergent for Hi 7^ and H2 7^ in asymptotic 
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solution J30j and ([31]) • For Hi = H 2 = 0, the radial flow 
velocities of both fluids remain finite at large x for n > 1. 
The radial profile of temperature in the hot ICM is T2 oc 
r 2-2/n £ Qr gjtjjgj. r _> q+ or r _> -(-oo. By mass conservation, 

the enclosed mass is continuous across shocks. 



5 DISCUSSION 

There are several relevant aspects that may be further taken 
into account in our theoretical model development. First, we 
have ignored the magnetic field permeated in the hot ICM 
which has been observationally inferred to range from a few 
fiG to several tens of fiG in the central region (e.g., Fabian 
1994; Carilli & Taylor 2002; Hu & Lou 2004). Cluster mag- 
netic field can affect the gas dynamics in a nontrivial man- 
ner and can provide valuable diagnostic information such as 
synchrotron emission and modified Sunyaev-Zel'dovich ef- 
fect (e.g., Hu & Lou 2004; Bagchi et al. 2006). To better 
understand behaviours of a magnetized hot gas and dark 
matter coupled by gravity, a magnetohydrodynamic (MHD) 
approach may be adopted (e.g., Lou & Wang 2007; Wang 

6 Lou 2007). Secondly, we have assumed the flow system 
to be spherically symmetric, which is a highly idealized sit- 
uation. Various instabilities may arise in dynamic flows to 
destroy spherical symmetry (Lou & Bai, 2007 in prepara- 
tion). For systems of merging galaxies or galaxy clusters, 
spherical symmetry is only a gross approximation. Thirdly, 
our current model requires that the mass density profile has 
the same form of power law scalings either at small or large 
radii. If the scaling parameter n can be adjusted in some 
proper way across a shock, this model may be adapted to 
simulate more diverse kinds of galaxy clusters. Finally, it is 
possible to take a distribution function approach for dark 
matter halo and fluid description for the hot ICM to model 
a galaxy cluster as a further development. Admittedly, this 
approach could be mathematically challenging. 

Compared with previous results, we note a few points 
regarding our model. First, our dynamic model of galaxy 
cluster wind differs from that of Yahil & Ostriker (1973). 
They considered a steady-state wind while the galaxy cluster 
wind of our model is dynamic and self-similar. Their gravi- 
tational potential is static and they did not model shock fea- 
tures occurred in galaxy clusters. Secondly, the case of n = 1 
for a single fluid sphere was discussed in Suto & Silk (1988). 
In addition to the major difference of two fluids coupled by 
gravity, our model with n = 1 differs from that of Suto & 
Silk (1988) in that we have specific entropy conservation 
along streamlines. In Lou (2005), the model of two isother- 
mal fluids coupled by gravity is discussed. Here, we treat 
two polytropic fluids coupled by self-gravity in more general 
situations. In particular, the quasi-static solution does not 
exist in an isothermal fluid (Lou & Wang 2006) . By choosing 
different values of parameters n, ji, and n and shock posi- 
tions, our polytropic model can be adapted to various astro- 
physical systems, including clusters of galaxies and globular 
clusters on much smaller scales. 

Due to possible galaxy cluster winds, baryons in galaxy 
clusters can flow out of galaxy clusters and stay in their pe- 
riphery. Meanwhile, dark matter halo contracts under grav- 
ity. As time goes on, such baryon gas cools down and may 
not be detectable in X-ray bands. This may naturally ex- 



plain the problem of 'missing baryons' for clusters of galax- 
ies. With this scenario in mind, we should develop obser- 
vational diagnostics to look for signatures of these missing 
baryons around clusters of galaxies on much larger spatial 
scales ( > Mpcs). In our perspective, the giant radio arcs 
recently discovered at 1.4GHz in galaxy cluster Abell 3376 
(Bagchi et al. 2006) and the strikingly similar radio arcs in 
galaxy cluster 3667 (e.g., Rottgering et al. 1997; Roettiger et 
al. 1999) are most likely large-scale magnetohydrodynamic 
(MHD) shocks and imply large-scale flows around clusters 
of galaxies (Lou & Jiang 2008 in preparation). 

As mentioned earlier, AGN activities are directly ob- 
served through intense electromagnetic radiation from nor- 
mal matter. For AGN activities in the central region of a 
galaxy cluster, accretion of dark matter onto a SMBH may 
be involved (e.g., Hu et al. 2006) although dark matter activ- 
ities cannot be directly detected. In this scenario, both ICM 
and dark matter coupled by gravity are active components 
of AGN activities and violent relaxation occurs in highly dis- 
turbed dark matter halo (Lynden-Bell 1967). In this sense, 
a dark matter halo gains gravitational energy and 'thermal' 
energy with higher velocity dispersions. Furthermore, shocks 
can also form in association with cluster merging processes 
also observed to occur in clusters of galaxies (e.g., the galaxy 
cluster IE 0657-56; Clowe et al. 2006). In such merging pro- 
cesses, dark matter is dragged along and the energy may be 
transferred to the dark matter. This mechanism has already 
been discussed by some authors (e.g., Shchekinov & Vasiliev 
2006; Knebe et al. 2002). 

Very recently, caustics in dark matter halo have been 
discussed by many authors (e.g., Natarajan & Sikivie 2007; 
Onemli & Sikivie 2007; Mohayaee et al. 2007). Caustics are 
consequence of collisionless cold dark matter in galactic ha- 
los or clusters of galaxies. Caustics are regions of infinite 
density in the limit that the DM particles have zero velocity 
dispersion (e.g., Natarajan & Sikivie 2007). Gravitational 
lensing effects have been proposed to detect the existence 
of such caustics (e.g., Onemli & Sikivie 2007). Compared 
with sharp drops of DM density in caustics, shocks in dark 
matter halo of our model are different in that the disconti- 
nuity of DM density across shocks is less drastic. In general, 
the density on the downstream side of shocks are ~ 2 — 10 
times of the density on the upstream side of shocks. How- 
ever, DM shocks share certain common features with DM 
caustics. They are both discontinuous surfaces in dark mat- 
ter halo and both may be detected through gravitational 
lensing effects at least in principle. 
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APPENDIX A: SELF-SIMILAR NONLINEAR 
ORDINARY DIFFERENTIAL EQUATIONS 



On the basis of reduced self-similar equations (|13p — (|16[) 
in the section of Model Formulation, the reduced radial flow 
speeds vi (x) and v 2 (x) and the reduced mass densities ai (x) 
and a 2 (x) are determined by four coupled first-order non- 
linear ODEs shown below. 



dai(x)/dx — A\(x)/T>i(x) , 
dvi(x)/dx — Vi(x)/T>i(x) , 
dct2(x)/dx — A2 (x)/T>2 (x) , 



(Al) 
(A2) 
(A3) 



APPENDIX B: THE QUARTIC EQUATION 

FOR INDEX PARAMETER K 



In equation (|26[) . the index k for the quasi-static solution of 
two gravity coupled fluids are determined by the following 
quartic equation 



Cfc.i k 4 + C* fe ,2 fc 3 + C fc , 3 k 2 + C kA k + C k , 5 = , (Bl) 



where the five coefficients Ck,i, Ck,2, Ck,3, Ck,4, and Ck,s 
are explicitly determined by the following expressions 



dv 2 (x)/dx = V 2 (x)/V 2 (x) . (A4) 

The explicit expressions of the two denominators T>\ (x) and 
T> 2 (x) are given by 

Vi(x) = (nx - vi) 2 - 11 a q 1 1+J1 - 1 x 2c ' 1 (nx - , (A5) 

and 

T>2(x) = (kiii — V2) 2 /k 

- K 2q2 - 1 ~f2a q 2 2+ ~ t2 - 1 x 2q2 {nnx ~ v 2 )" 2 . (A6) 

The explicit expressions for the four numerators Ai(x), 
Vi(x), A2 {x) and V2{x) are given by 

Ai{x) = a 1 L 1 a q 1 1+ ' 11 - 1 x 2qi (nx - Wi) 9l-1 (3n - 2) 



(A7) 



Oil , \ . OL2 I V 2 

+ -, r ( nx — Vi 1 + -, r nx 

(3n-2) v ; (3n-2)V K 



+ (n- l)vi 



2(x — vi)(nx — v\) 



Vi(x) = qia qi+11 - 1 x 2qi (nx-v 1 ) qi {?,n-2) 



2 7l Q? 



(nx — Vi) qi (x — v\) 



ai(nx — vi) 2 a 2 (nx — vi)( v 2 
+ (3n - 2) + (3n - 2) \ nX ~ ~ 
+ (n — l)(nx - Vi)vi ; 



(A8) 



A 2 (x) = a 2 { q2K q2 a q 2 2+ ~ ,2 ~x' q2 (Knx - v 2 ) q2 ~ (3n - 2) 



a 2 (nnx — 1)2) KO>i(nx — Vi) 
(3n - 2) ' (3n - 2) 

2(kx — v 2 )(Knx — v 2 ) 
+ (n — 1)V2 



(A9) 



V2(x) = q 2 K 



2 12 „,92+72~l 2<2 2 



■ a 



x zq2 (nnx - v 2 f 2 (3n - 2) 



-2n 2q2 1 72«2 2+72 1 x 2q2 1 (nnx — v 2 ) q2 (kx — v 2 ) 
a 2 (nnx — v 2 ) 2 nai(nx — vi) 



(3n - 2) + ' (3n - 2) 
+ (n — l)(nnx — v 2 )v 2 . 



-(nnx — v 2 ) 



(A10) 



Here, all relevant dimensionless parameters are defined in 
the main text. With proper asymptotic solutions at large 
and small x and shock conditions across the sonic critical 
curves, we can construct global semi-complete self-similar 
solutions with shocks to model large-scale dynamics in a 
galaxy cluster involving hot ICM and dark matter halo. 



(B2) 



< 1 1 = f 1." - 8)n 3 7l72 ( 1 • x) ' T" 



(B3) 



C fc ,3 = I - 3n 4 7i 7 2 + ra 2 (4 - 2n)(q 2 + 72 - 1) 

x [(2-3n) 7 i-(4-2n)(gi-l)] 

+ n 2 (4 - 2n)( qi + 71 - 1)[(2 - 3n) 72 

- (4 - 2n)(q 2 - 1)] + (4 - 2n) 2 n 2 ( Ql - l)(q 2 - I) 

+ (4-2n)V(gi +71 - l)(q 2 + 72 - f)| 



+ (3n-2)(4-2n)n 2 7i( f + 



+ (3n-2)(4-2n)n 72 1 + 



A. 
A, 

M 
A2 



(B4) 



C fcl4 = 2 1 [ 72 n + (2 - n)(q 2 - l)]n 3 7 i 

+ [ 7 m + (2 - n)(qi - f)]n 3 7 2 
- 4[(2 - n)(qi - 1) + 7 in][ 72 n + (2 - n)(q 2 - l)]n 2 
+ 2n(q 2 + 72 - 1)(2 - n)[n 7 i + (2 - n) 2 ( qi - I)] 
+ 2n( Ql + 71 - 1)(2 - n)[nj 2 + (2 - nf(q 2 - I)] 

+ 8(2 - nf(n - l)n(q 2 + 72 - l)(gi + 71 - 1) } 



+ 2(3n - 2)(3n - 4) (2 - n)n7r f + 



M 
Ai 

+ 2(3n - 2)(3n - 4) (2 - n)n 72 ( 1 + 4^ 

jT-2 



(B5) 
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r. , = 1C!„ - 2)12 - n)[ir - (in + ( 1 + ^ 
+ 4(3n - 2) (2 - n)(n 2 - 6n + 4) 72 f 1 + 4 1 

V A 2 

+ 4[ 7l n + (2 - - l)][ 72 n + (2 - 7i)(g 2 - 1) 

+ 8[ 7 m + (2 - n)(qi - 1)](2 - n)(g 2 +72 - 1) 

xn(n-l)(l + |)(l + f 
+ 8[ 72 n + (2 - 7i)(g 2 - 1)](2 - n)(gi + 7 i - 1) 

+ 16(n-l) 2 (2-n) 2 (g 2 + 72 -l) 

K fe+ *-l>(l + *)(l + * 

-4(3n-2) 2 (2-n) 2 (g 1 -l)('l + ^ 2 - 



4(3n-2) 2 (2-n) 2 (g 2 -l) 1 + 



A 2 



(B6) 



where the two coefficients Ai and A 2 in the static SPS so- 
lution (|17|l are determined by equation (|18|l . In short, once 
the four parameters (n, 7 i, 72 , k) are specified, we should 
be able to determine the static SPS solution (|17p and the 
corresponding quasi-static solution (f2"4")l — (|2"T|) . 



APPENDIX C: QUASI-STATIC SOLUTIONS 
WITH COMPLEX K 

For a complex k root with k = k± + ifc 2 with ki and fc 2 being 
the real and imaginary parts of k (note that the value of k 
is the same for both fluids coupled by gravity), the four pa- 
rameters Li, L 2 , Ni and iV 2 (the coefficients for radial flow 
speeds and mass densities at small x of quasi-static solutions 
defined in equation (|24p and Q25p) all become complex with 
Li = Li,i +iLi, 2 , L 2 = L 2 ,i + iL 2 ,2, Ni = 2Vi,i +iA r i, 2 , and 
N2 = Ns,i + iA^s^ 1 where the real and imaginary parts are 
all explicitly written out. As the static SPS solution (|17p is 
real, we only take the real part as the quasi-static solution, 
namely 

Vi(x) = Re(LiX k ) = Re[LiX kl exp(ifc 2 lna;)] 

= x kl [Li,i cos(fc 2 lna;) — Li, 2 sin(fe 2 lna;)] , (CI) 



cti(x) 



Re\ Aix~ 2/n +NiX 



k-X-2/n 



Re 



' 2/n + x - 2/n - 1+kl Ni exp(ife lna;) 



^a;- 2/ " + a;- 2/ "- 1+fcl 



x [iVi.i cos(fc 2 In x) - Ni,2 sin(fc 2 lna;)] , (C2) 



where 



n 2 [(fci - l) 2 + fc§]JVi, a = (fci - l)[2(n - 1) 

+nk 1 ]L 1>1 A 1 + (3n - 2)fc 2 Li, 2 Ai + nfe|Li,iAi , (C3) 

Ni,s = Ni,i(ki-l)/k2 + L 1>2 Ai/n 



7i 2 [(fci - l) 2 + fc|]«JV 2 ,i = (&i - l)[2(ra - 1) 

+ nk 1 ]L 2 ,iA 2 + (3n - 2)fc 2 L 2 , 2 A 2 + nfc|L 2 , 1A2 , (C5) 

A*2,2 = N 2 ,i(ki - l)/ft 2 + L 2 , 2 A 2 /(ft;n) 

-[2(n - 1) + nfci]L 2 ,i^ 2 /(Kfc2n 2 ) . (C6) 

As equation (|26p is a quartic equation for k involving real 
coefficients, it has either real roots or pairs of complex con- 
jugate roots. With free parameters Li,i and Li, 2 , we can 
always choose the imaginary part of a complex roo10 k to 
be k 2 > 0. 

For this type of quasi-static solutions with complex k, 
there are asymptotic oscillations as x — * + (see also Lou 
& Wang 2006 for such oscillatory solution behaviours in a 
single fluid). 



APPENDIX D: ASYMPTOTIC SERIES 
SOLUTIONS AT LARGE X 



As stated in subsection [231 f° r a given value of n, the asymp- 
totic series solution of equations (|13p — (|16p at large x is in 
the form of 



a 1 (a;)^L 1 ar 2/n +/ 1 ar 3/n + ... , 
a 2 (aO^L 2 a;- 2/n +/2a;- 3/n + ... , 

( \ rr — 1 / n+1 , — 2/ n+1 . 

vi (x) — ► Hix ' +Gix 1 +. 

r \ tt — 1 / n+1 , r~i —2/ n+1 . 

V2(x) ^ H2X ' T + G 2 a; ' +. 



(Dl) 
(D2) 
(D3) 
(D4) 



The four coefficients 7i, I2, G\ and G2 can be expressed in 
terms of Ei, E2, Hi and H 2 as 



Gi = (1 Tl) Hl + 2(2 - n)El 1+11+1 n qi - 1 
n 

(E 1 +E 2 ), 



G 2 = 



(3n - 2) 



(1 - n) g 2 2 
n k 
nn 

~ (3n - 2) 



+ 2(2 - n^^Kp+^+V 2 " 1 



(Si + E 2 ) 



h = 



3(1 - n) H2E2 



n k 



(D5) 

(D6) 
(D7) 
(D8) 



-[2(n - 1) + nk 1 ]L hl A l /(k 2 n) 



(C4) 



For solutions with n > 1 and Hi = H2 ~ 0, the asymptotic 
series solutions finite at large x become 

Qi(a;)^£ia;- 2/n +F 1 a;- 4/ ' l+1 + ... , (D9) 



6 Parameter k2 is the imaginary part of a complex k index. If 
k 2 < 0, for the k 2 that appears in function cos(x), it does not 
matter to choose it to be — k 2 . For k 2 that appears in function 
sin(a;), as long as Li j2 being a free parameter, we can always 
choose it to be — L^2 and it follows that Ni 2 becomes —Ni }2 . In 
this manner, solution for — k2 remains the same as that for 
In conclusion, there is no loss of generality to choose k 2 > to 
represent all possible solutions with free parameters Li. 
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ai2 (a) -» E 2 x- 2/n + F 2 x~ 4/n+1 + ... , 
vt(x) -> GiaT 2/n+1 + + . . . , 

v 2 (x) - G 2 z- 2/n+1 + J D 2 ^- 4/n+2 + . . . . 



(D10) 
(Dll) 
(D12) 



The six coefficients Gi, G 2 , Fi, F2, D\, and D2 are deter- 
mined by specified values of E\ and E2 below 



G 1 = 2(2 - r^Ff+^+V 1 " 1 - n(i?1+Jg2) 



(3n - 2) ' 



G 2 = 2(2 - n)« J 



392-l^j92+72 + l„<32 



1 92 _1 /tn(-Bl + £2) 



(3n - 2) 



Fi = (4 - 3n)(2 - n)Sf +TI n ?1 - 2 
(3n - 4) 



2(3n - 2) 



E l (E 1 +E 2 ) 



F 2 = k 3?2_2 (4 - 3n)(2 - n)E 2 
+ (3n-4) 



2(3n - 2) 



F 2 (Fi +F 2 ) 



(D13) 



(D14) 



(D15) 



(D16) 



£>i = 



(3-n) 



(2-n)(gi +71 - l)S 1 1+71 - 2 n 91 - 1 F 1 



n(Fi + Fi) 



D 2 = 



(3-n) 



2(3n - 2) 



(2-n) K 392 - 1 ( (?2 + 72 -l) 



(D17) 



x e: 



92+72-2„92 



1 „ «n(Fi + F 2 ) 



2(3n - 2) 



(D18) 



L'Hospital rule to determine the eigenvalues of «l(a;) from 
equation (fT4|) 



d«i(a;) Vj(a;) 
da; _ £>i(ai) 



(E4) 



This equation appears to be a quadratic algebraic equation 
in terms of v[(x), namely 



Gv,i Wi(x)] 2 + a,a v[(x) + a, 3 = 



(E5) 



where the three coefficients C V) i, C Vt 2, C Vl z are explicitly 
defined by 

CVi = 1 + a? +11 ~ 1 1 2 1 x 2qi {nx - Vi)" 1 - 2 , (E6) 

C v , 2 = a^+^^j^-^nx ~ v 1 ) qi - 2 {4(j 1 - l)vi 

+[3n(«i + 1) - 2(9i + 71 + 1)M - 1 , (E7) 



C v ,3 = 2 a f + ^- L x^"{nx^v l ) qi - A ^ n {2 11 -l)vt 

+ [n(3gi7i - 2) - 2[ 7 i (qi + 71 + 1) - 2]] xvi 

- (n-2)[n-3n 9 i + 2(,3i+7i-f)].T 2 
na.2 Qi (n — 2 + 2i>i /a;) 



+ 



(3n - 2) (3n - 2) 



£¥2^2 



Q 2 W 2 nxQ 2 



(3n - 2)k (3n - 2)/s (3n - 2) 



(E8) 



We can then readily solve for the two eigenvalues of v[(x), 
corresponding two eigendirections. The two possible eigendi- 
rections across the corresponding sonic critical curve of fluid 
2 can be computed in the same manner. 



APPENDIX E: EIGENDIRECTIONS ACROSS 
THE SONIC CRITICAL CURVE 

For a specified pair of (a 2 , V2) at a given x, we now determine 
the eigendirections across the sonic critical curve of fluid 1. 
The first equation in nonlinear ODEs (Equation (I13p ) gives 



dai 
dx 



a 1 



dx 



n (x- vi) 

I Ql 



(nx — Vi) 



(El) 



We take n > 2/3 and thus nx — v > in this paper. From 
'Di(x) — 0, we obtain 



vi — nx 



7 ix 291 af + 71 



"I 



1/(2-91) 



(E2) 



From Ai(x) — 0, the value of «i on the sonic critical curve 
at x is determined by the following equation 



01 



3?i-2 



+ n- 



2n-4 



2(7i"? 

a 2 
' (3n - 2) 



71 



91+71-1 



(7i«? 



91+71-1 2giNl/(2- 8l ) 



2/(2-91) 



,(59i-2)/(2-q!) 



t>2 
K 



n(n — l)x 



(E3) 



In general, fluid 2 is not on the sonic critical curve at this 
position x. Therefore, a' 2 is determined by equation (|15p 
while v' 2 is determined by equation (|16[) . We now apply the 



APPENDIX F: SHOCK CONDITIONS 



As stated in subsection of self-similar shocks [2^61 all variables 
with subscript {u, i} refer to variables on the upstream side 
of a shock for fluid i, while all variables with subscript {d, 1} 
refer to variables on the downstream side of a shock for fluid 
i. Convenient variables V U: i and Td,i are defined by 



r«t, 



Vd,i 
Xd,i 



Parameter A; defined by 



Xi = 



1/2 



(Fl) 



(F2) 



(F3) 



is essentially the sound speed ratio of downstream and up- 
stream sides of a shock and should be greater than 1 for 
entropy increase on the downstream side. According to the 
dimensionless equation of state (see equation (HJ), the di- 
mensional equation of state is 

Pi = ^ I (1 - 39l/2) (47r) 7 '- 1 G 7 '- 1+9 '(3n- 2) 9 *M«*p 7 ' . (F4) 

Then shock equations (equation of mass conversation, radial 
momentum conservation equation and energy conservation 
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equation ([52]) — J54J)) can be rearranged into dimensionless 
self-similar forms of 



otu,i{nx Uti — v u a) = ctd,i{nx U:i — XiV d<i ) 
= ct dA \i(nx d>i — v d>i ) , 

\ a 'd,t' 1i x d 9 J ( nx d,i - v d ,i) 9i + ad,i(nxd,i - «d,i) 2 ] 



^uA 



(F5) 



(F6) 



reduced enclosed mass, the polytropic sound speed in fluid 
i can be written as 



(G3) 



Since {KyjK^) = X2/X1 = k, the sound speed ratio of 
fluid 1 to fluid 2 is given by 



ai 

«2 



72 



.1-92 111 -91-93 [ - UXl Vl > 



„(92+72-l)/2 



(nx 2 - V 2 ) q2/2 



(G4) 



9i+7i-l 



+ (7< _i) g.i.H™*. 



2 ' (7i 

(nx Uti - v uA ) 2 /2 

+7i a u,i Ji ^ 1 ^K'nxu,i - v u ,i) qi /{H - 1) 



(F7) 



With Ti and T2 for the two coupled fluids, we can rewrite 
the above three equations in more symmetric forms, namely 



v 9i+7i„3tJi-2 r?i 



^dli + ad,iT dii 



_ n ,«i+7i„3qr i -2 r , 4 



' OCu,iJ^u A 
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2 ( 7l -l) d > 1 d ' 1 



2' + ( 7 



7? n/ 9i+7i-l r 3qi-2-pg i 



(F8) 



(F9) 



(F10) 



From above three equations (|F8|I — (|F10|) . F u a upstream of a 
shock can be determined by the variables on the downstream 
side of a shock from the following quadratic equation 



(7* + 1) p2 



2 7 i 



<K+7i-l r <7i-l 3 9l -2 p 
u d,i 1 d,i x d,i Tl-ii,! Ilu.! 



+Q ,9i+7i-l r 8» 3 9j -2 (7 



27« 



l) r 2 
— 1 d,i 







(Fll) 



where subscript i = 1, 2 correspond to fluid 1 and 2, re- 
spectively. One root of quadratic equation (|F1 1 f) is a trivial 
solution T u ,i = Fd,i, which is omitted. The other root of 
physical relevance is given in the subsection of self-similar 
shocks 12.61 



The sound speed varies with radius and the sound speed ra- 
tio at different radii can be computed according to equation 

(Gl- 
In the regime of x — * + , we take the quasi-static solu- 
tion of our model and the asymptotic behaviour of density is 
determined by equation (|25[) (velocities are small compared 
with x). Thus when x — * + , the sound speed ratio has the 
asymptotic behaviour 



ai V 72 / 



71 Y' * K l-3q 2 /2 Ay 



A 



(92+72-l)/2 



(G5) 



As the coefficients of static polytropic solution Ay and A2 
are determined by equation (|18|l . the sound speed ratio at 
small radius is 



01 _ /71 n q2 \ i/ 2 
0,2 V 72 n qi . 



(G6) 



Therefore in the limit of x ~ > + and for the quasi-static 
solution, the sound speed ratio is independent of n and is 
only dependent on n and on the two polytropic indices 7,. 

In the limit of x —* +00, the asymptotic behaviour of 
density is determined by equation (128[l and equation (|29[1 
(velocities are small compared with x), so the asymptotic 
behaviour of sound speed ratio at large x is 



01 
a 12 



7iy/\i- 3 « 2 /2 E± 

~12> 



(91+71-l)/2 



(92+72-l)/2 



E. 



(G7) 



where Ey and E2 are two coefficients determined by equa- 
tions ([28]) and Ip9]l. 

Two specific examples for the radial profile of polytropic 
sound speed ratio a 2 /ay are shown in Fig. [5] and Fig. 1151 in 
the main text. 
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